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ABSTRACT 


This research contributes to multiple spacecraft control by developing an 
autonomous distributed control algorithm for close proximity operations of multiple 
spacecraft systems, including rendezvous and docking scenarios. The proposed control 
algorithm combines the efficiency of the Linear Quadratic Regulator (LQR) and the 
robust collision avoidance capability of the Artificial Potential Function (APF) method. 
The LQR control effort serves as the attractive force toward goal positions, while the 
APF-based repulsive functions provide collision avoidance for both fixed and moving 
obstacles. The combination of the LQR and APF control logics, referred to as the 
LQR/APF control algorithm, yielded promising results as demonstrated by the numerous 
multiple spacecraft maneuver simulations reported in this dissertation. 

In order to validate the proposed control approach, a multiple spacecraft model 
validation and visualization technique was developed using a versatile MATLAB- 
Satellite Toll Kit (STK) interface to propagate the spacecraft models, compare against 
STK generated ephemeris, and animate for analysis. The MATLAB-STK interface 
efficacy was demonstrated during the evaluation and analysis of the innovative LQR/APF 
multiple spacecraft control algorithm. 

The LQR/APF multiple spacecraft close proximity control algorithm was 
developed, refined, and thoroughly simulated using high fidelity six Degree of Freedom 
(DOF) spacecraft models. In order to evaluate the stability and robustness of the control 
approach a Monte-Carlo simulations set was run. The LQR/APF control algorithm was 
further evaluated by virtual hardware-in-the-loop implementation at the NPS Spacecraft 
Robotics Laboratory. The laboratory hosts the Autonomous Docking and Spacecraft 
Servicing testbed which allows for on-the-ground testing of close proximity multiple 
spacecraft control concepts. 
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EXECUTIVE SUMMARY 


As spacecraft technology has advanced, simultaneous control of multiple 
cooperative spacecraft has become a desired mission capability. The first generation of 
spacecraft were individually stabilized in orbit by manual control from ground stations. 
The second generation of spacecraft were placed in stable constellations by ground 
stations, with some basic controls automated on-board the spacecrafts. The next 
generation will be required to reliably perform autonomous multiple spacecraft close 
proximity operations while avoiding collisions. Therefore, a need exists for robust and 
efficient automated and distributed control of multiple spacecraft for emerging servicing 
missions, involving simultaneous rendezvous and docking scenarios. 

In this research a novel multiple spacecraft close proximity control algorithm was 
developed, refined, and thoroughly simulated using high fidelity six degree of freedom 
spacecraft models. The developed control algorithm combines the dynamic optimization 
of a Linear Quadratic Regulator (LQR) and collision avoidance capability of Artificial 
Potential Field (APF) approaches. Development and evaluation of the LQR/APF control 
algorithm was supported by the realization of a MATLAB and Satellite Tool Kit (STK) 
simulation interface technique. This versatile MATLAB-STK simulation interface 
allowed for both multiple spacecraft model validation and simulation visualization. The 
LQR/APF multiple spacecraft control algorithm was further evaluated by virtual 
hardware-in-the-loop (VHfL) configuration at the NPS Spacecraft Robotics Laboratory 
(SRL). The VHfL structure utilizes independent processors which simulate spacecraft 
and interact as multiple spacecraft. 

Analysis of numerous close proximity maneuvers proved the LQR/APF to be both 
effective and efficient in conducting simultaneous spacecraft missions. The LQR/APF 
avoids actuator saturation while avoiding both stationary and moving obstacles. Monte 
Carlo simulations showed the multiple spacecraft control to be both stable and robust, ft 
also established a parametric baseline for future multiple spacecraft close proximity 
control algorithms, by evaluating collision avoidance requirement, fuel efficiency, and 
maneuver duration objectives. The LQR/APF control algorithm’s desirable performance 



gives spacecraft designers and mission planners a useful means of developing and 
forecasting maneuvers. Finally, successful implementation of the multiple spacecraft 
control algorithm in a VHIL configuration paves the way for future terrestrial and orbital 
hardware testing. 
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I. INTRODUCTION 

A. MOTIVATION AND BACKGROUND 

As spacecraft technology has advanced, simultaneous control of multiple 
cooperative spacecraft has become a desired mission capability. Therefore, there is a 
need to develop an autonomous distributed control algorithm for multiple spacecraft 
during close proximity operations. The first generation of spacecraft were individually 
stabilized in orbit by manual control from ground stations. The second generation of 
spacecraft were placed in stable constellations, and formations by ground stations, with 
some basic controls automated on-board the spacecrafts. The next generation of 
spacecraft will be required to reliably perform autonomous close proximity operations, 
including multiple spacecraft dispersion, rendezvous, and docking maneuvers. 

There are numerous mission concepts that involve the convergence of multiple 
spacecraft in close proximity. At present, these missions are parameterized through pre¬ 
determined (a priori) orientations and trajectories, and executed with centralized manual 
control. This approach is extremely cumbersome with respect to time and computational 
expense, relies on high communication capacity between spacecraft, and allows for no 
dynamic reconfiguration in spacecraft control. While, large spacecraft formation tracking 
and keeping has received a great deal of study, research in the area of practical multiple 
spacecraft close proximity operations is limited [1]. Current spacecraft rendezvous and 
docking require that the spacecraft cooperate and are deliberately designed to work 
together. Approaching spacecraft usually use tracking based algorithms, with advanced 
sensors and processors, to approach the desired position. Additionally, these close 
proximity spacecraft maneuvers are typically manually operated. The need for advanced 
sensors, processors, and actuators limits the application of previous spacecraft control 
algorithms to a wide variety of spacecraft platforms and mission. 

These issues could be addressed with a computationally efficient, robust, 
distributed control algorithm allowing for multiple spacecraft close proximity operations, 
with the capacity for dynamic reconfiguration for collision avoidance. Research and 
experience with terrestrial based robots have matured the application of iterative artificial 
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potential field (APF) based control algorithms [2][3][4][5]. These APF algorithms are 
geometrically based with converging minimization searches along the direct Steepest 
Descent path. The APF have been refined to include obstacle repulsion potentials which 
allow for robots to avoid collisions in moderately dynamic environments [6] [7]. 
However, these algorithms have not been widely applied to the unique challenges 
presented in spacecraft control. The simplicity of the potential field based control 
algorithms are a good match for spacecraft application with limited proximity sensors and 
processing capability. Previously proposed spacecraft potential field based controllers 
have been task and platform specific and not robust in the full range of possible close 
proximity operations. Also, fuel efficiency and optimization has been limited to 
maintaining spacecraft formations. The consideration of efficiency while maintaining 
collision avoidance in close proximity operations has been limited, and typically requires 
dramatic increases in computation or centralization. 

Efficiency concerns have driven spacecraft control system designers toward 
optimal control algorithms [8]. However, the restricted computation capability onboard 
the current generation of spacecraft forces trade-offs between such algorithms and simple 
feedback architectures which can be hosted locally. An original solution to the efficiency 
demands and the collision avoidance capability required for the emerging autonomous 
close proximity operations is the combination of optimal control and geometric collision 
avoidance. An iterative Linear Quadratic Regulator (LQR) control algorithm can be used 
as driving force toward the desired goal, replacing the more conventional force due to the 
attractive potential. The LQR has the advantage of incorporating the spacecraft 
linearized relative dynamics by utilizing variable state and control effort gain matrices in 
order to solve for an optimal cost solution at each iteration. These allow variation for 
state weighting during convergence toward the goal. 

The fusion of a LQR, including relative dynamics, and an APL based collision 

avoidance capability yields a new and promising multiple spacecraft control algorithm. 

The LQR control effort serves as the attractive force toward goal positions, while the 

APF-type repulsive functions provide collision avoidance for both fixed and moving 

obstacles. This LQR/APL multiple spacecraft proximity control algorithm offers 

desirable performance in a robust close-proximity, and establishes a baseline for fuel 
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efficiency while maintaining collision free operations. Critical evaluation of the multiple 
spacecraft control algorithm utilized high fidelity six DOF spacecraft models, allows 
assessment with realistic spacecraft dynamics and constraints. The multiple spacecraft 
close proximity control algorithm was evaluated for several multiple spacecraft emerging 
maneuvers, which may require gathering, rendezvous, and docking operations. 

This dissertation includes the requirement, development, simulation, refinement, 
evaluation, and analysis of an autonomous distributed control algorithm for multiple 
spacecraft during close proximity operations. Literature reference and review material is 
presented in each applicable dissertation chapter. This ensures topical material is 
presented as research concepts are presented and developed. Chapter I closes with a 
summary of the research and its contributions. The mission profile and requirements are 
presented in Chapter II. The space environment and a high fidelity six DOF spacecraft 
model are developed in Chapter III. The equations of relative motion between multiple 
spacecraft in close proximity and the typical rendezvous control are introduced in 
Chapter IV. The close proximity spacecraft control algorithm is developed in Chapters 
V-VII. A detailed evaluation of the LQR/APF performances during convergence, rally, 
rendezvous, and docking maneuvers is presented in Chapter VIII, followed by a Monte- 
Carlo method analysis of the heuristic control algorithm in Chapter IX. Finally in 
Chapter X, a virtual hardware-in-the-loop (VHIL) implementation of the multiple 
spacecraft control algorithm is discussed. 

B. RESEARCH GOALS 

This research is intended to advance the field of multiple spacecraft control by 
developing an autonomous distributed control algorithm for multiple spacecraft in close 
proximity operations, including simultaneous rendezvous and docking. The architecture 
of the controller is designed to convolve the collision avoidance capacity of APF 
approaches with the LQR capacity to address dynamic platform constraints. The 
developed approach, referred to as the LQR/APF control algorithm, is an iterative 
feedback based algorithm which allows for efficient and timely proximity spacecraft 
maneuvers. The LQR/APF control algorithm combines efficient LQR performance. 
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utilizing linearized relative dynamies, and APF with geometric collision avoidance. This 
convergence of efficiency and collision avoidance into a control algorithm is an enabler 
for future spacecraft missions. 

The developed control algorithm was thoroughly simulated with a high fidelity 
six degree-of-freedom (DOF) spacecraft model operating in 3D space. Next, the 
simulation model was fully developed with consideration of both perturbation forces and 
torques. The control algorithm takes into account considerations on realistic actuators 
and sensors performances. Finally, the control algorithm was tested and evaluate by 
VHIL implementation in the dedicated Spacecraft Servicing and Robotics Laboratory 
located at the Naval Postgraduate School (NPS). This allows for validation of the control 
algorithm with independent processors simulating spacecraft performance based on 
limited and incremental state information. The laboratory also hosts the Autonomous 
Docking and Spacecraft Servicing test-bed, whose task is to simulate on-the-ground the 
navigation and control of the docking between multiple free-flying small spacecraft. 


C. MAIN CONTRIBUTIONS 

Specific research contributions include: 

• Development of a MATLAB-Simulink interface with STK which allows 
for multiple spacecraft model validation and dynamic environment 
visualization. This developer friendly tool proved critical for engineering 
analysis of control algorithm performance. The use of STK, which is a 
key spacecraft industry standard, allows for clear presentation of 
developed control algorithm performance to the spacecraft field. 

• Development of a multi-purpose APF based control algorithm which 
performs close proximity operations in the spacecraft environment with 
realistic spacecraft constraints. This algorithm was refined with the 
addition of collision avoidance capability to be effective in a broad range 
of spacecraft maneuvers involving multiple spacecraft and obstacles. 

• Development of an iterative LQR control algorithm, with variable gain 
matrices, which performs close proximity operations in the spacecraft 
environment with realistic spacecraft constraints. 

• Combination of the LQR and APF control concepts. The iterative LQR 
generated control effort serves as the attractive force toward goal 
positions, while the APF repulsive functions provide collision avoidance 
for both fixed and moving obstacles. This algorithm was further refined 
for a full range of maneuvers involving stationary and moving obstacles. 
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• Thorough comparison of maneuver performanee, effieieney, and duration. 
This researeh established a parametrie baseline for future multiple 
spaeeeraft elose proximity eontrol algorithms, by evaluating eollision 
avoidanee requirement, fuel effieieney, and maneuver duration objeetives. 

• Implementation and evaluation of the eontrol algorithm in a VHIL 
eonfiguration. Utilizes independent proeessors whieh simulate spaeeeraft 
and interaet as multiple spaeeeraft. 
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II. OVERVIEW OE MULTIPLE SPACECRAET MISSIONS 


A. MULTIPLE SPACECRAFT CONTROL 

The controlled spatial interaction of systems involving multiple vehicles, robots, 
aircraft, or spacecraft in close proximity is a complex task. How multiple vehicles are 
controlled depends on a wide range of operating environments, missions, objectives, and 
organizational combinations. Developing systems consisting of multiple autonomous 
vehicles that cooperatively perform a task or behavior is of paramount importance to the 
engineering community [9]. The distinctive space environment makes for particularly 
interesting and challenging control algorithm requirements based on orbital dynamics, 
communication limitations, and tight spacecraft sensitive to disturbances and constraints. 
The dynamics of the orbital spacecraft is affected by disturbances and perturbations. 
Traditionally, individual spacecraft maneuvers have been based on fuel efficiency 
requirements, since fuel is a critical resource in the life of a spacecraft. During large 
early staging maneuvers, the conservation of fuel is generally more critical then the 
timeliness of any particular task. Once in the desired mission orbit, the only regular 
translational maneuvers which most individual spacecraft perform are for general station 
keeping. Also, due to the relatively large distances between typical orbital objects 
collision avoidance maneuvers are rarely necessary and seldom executed. However, due 
to spacecraft technology improvements, there is now a greater desire to control multiple 
spacecraft in close proximity operations. These multiple spacecraft, close proximity 
operations require collision avoidance while each spacecraft executes the desired close 
proximity task. 

Advances in technology continue to decrease the size, weight of components 
while increasing the capability of payloads, sensors and processors. This trend has 
enabled the space industry to decrease the size of spacecraft. Small satellites with mass 
less than a few hundred kg are becoming more common. Launch opportunities are more 
readily available for smaller spacecraft since more launch vehicle can support them. 
Many of these relatively small satellites can be deployed into orbit from the same launch 
vehicle. Once in orbit the spacecraft will need to disperse or converge, relative to each 

other depending on their mission. As spacecraft get closer to each other the execution of 
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collision avoidance, rendezvous, and doeking maneuvers over a short timeline may take 
priority over optimizing individual spaeeeraft's fuel eonsumption. During these elose 
proximity maneuvers, the desires are to maintain effieient fuel management while 
aeeomplishing the maneuver objeetive quiekly. The longer eaeh spaeeeraft stays in 
transition during elose proximity the more preeise station keeping is required; therefore 
more propellant may be used. 

B. MULTIPLE SPACECRAFT MISSIONS 

Control algorithms for multiple spaeeeraft need to take the similarity and 
differenees of eaeh spaeeeraft into eonsideration, wtiieii is diseussed in detail in Chapter 
II.B.l. In this researeh eooperative homogeneous spaeeeraft with similar sensors and 
eontrol algorithms are eonsidered. Additionally, it is useful to define the phase of the 
spaeeeraft mission in relationship to the distanees from other spaeeeraft or the goal 
position. A diseussion of relative orbital mission phases is presented in Chapter II.B.2. 
Also, the organizational grouping of spaeeeraft depends on the manner in whieh the 
spaeeeraft maintain position with respeet to eaeh other. The various organizational 
groupings of spaeeeraft are diseussed in Chapter II.B.3. Preeision eontrol may be desired 
in some applieations, sueh as pointing spaeeeraft with high resolution imagery payloads, 
and loose eontrol may be desired in other applieations, sueh as station keeping of 
eommunieation satellites. Mission motivations for autonomously eontrolled and 
synehronized multiple spaeeeraft in preeise spatial eonfigurations are numerous, and 
inelude interferometry, eommunieations, and power generation. These eooperative 
spaeeeraft may need to be able to eonverge or diverge in a safe and effieient fashion. For 
servieing missions, it may be neeessary for two or more spaeeeraft to rendezvous, or even 
doek. 

First, multiple spaeeeraft interferometry is based on the idea that the preeise 
eontrol of the relative position and orientation of multiple spaeeeraft payloads eould 
result in performanee equivalent to a mueh larger spaeeeraft payload. Multiple spaeeeraft 
are positioned in order to form a distributed array. Distributing spaeeeraft payloads (e.g., 
sensors) into preeise spatially eonfigurations, or assembling multiple spaeeeraft on orbit, 

require higher level eontrol algorithms. Required eontrol algorithm eapabilities inelude 
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the establishment and maintenanee of a given formation and any possible reeonfiguration 
of the formation. Formation maintenanee eontrol has been given its due of attention [10] 
[ 11 ], 

Seeond, reoonfigurable eooperative formation mission eoneepts are based on 
eolleetives of small spaeeeraft working in unison to perform a greater funetion. For the 
purpose of spaeeeraft formation eontrol it is generally assumed that eaeh spaeeeraft is 
relatively homogeneous and exeeutes eentralized eontrol strategies in order to minimize 
fuel eonsumption during formation flight. Eaeh small spaeeeraft eoordinates with others 
in the group and shares proeessing, eommunieation, and payload or mission funetions. 
Defense Advaneed Researeh Projeets Ageney (DARPA) is eurrently studying the idea of 
fraetionated spaeeeraft arehiteetures whieh deeompose the overall spaeeeraft funetion 
into a free-flying network of eomponent modules. The traditional monolithie spaeeeraft 
might be replaeed with a group of smaller spaeeeraft interaeting wirelessly [12]. 
Individual spaeeeraft may be required to perform elose proximity maneuvers 
autonomously while the mission or eonfiguration ehanges. These fraetionated, modular 
spaeeeraft have advantages over traditional monolithie spaeeeraft, as analyzed and 
assessed in [12] [13]. Upgrading the group ean be done iteratively in order to inerease 
overall performanee and mission duration. This ean substantially enhanees the 
flexibility, responsiveness, robustness, and lifeeyele of the overall spaeeeraft funetion. 

Third, rendezvous mission eoneepts are usually based on a eooperative spaeeeraft 

approaehing another spaeeeraft within a eommon spatial region. The goal position is 

usually oeeupied by a eooperative spaeeeraft in the same orbital plane as the 

maneuvering spaeeeraft. The Spaee Transportation System (STS), often referred to as the 

Spaee Shuttle, and the ISS are often shown booking in orbit. The Spaee Shuttle is 

astronaut oontrolled as it books to the ISS. This prooess is not autonomous, but there is 

autonomous booking of Russian Soyuz and Progress spaeeeraft with the ISS. The 

prooess of automated spaeeeraft booking was pioneered by the Soviet Union; however the 

automated system oooasionally fails to oomplete the task. Aooording to NASA, ourrent 

state of the art Russian automated rendezvous and booking systems have a ourrent failure 

rate of approximately 10-15 % [14]. As a result, ISS’s Zvezda module is equipped with 

the Russian built Telerobotioally Operated Rendezvous Unit (TORU) manual booking 
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system which can be operated by cosmonauts [15]. However, manual rendezvous control 
can not be relied upon since human space flight is very dangerous, impractical, and 
unnecessary for most on orbit missions. The European Space Agency (ESA) is currently 
developing the Automated Transfer Vehicle (ATV) to automatically dock with the ISS. 
However, close proximity operations of spacecraft are often complicated by momentary 
communication and autopilot navigation failures. 

The motivation to rendezvous is not limited to manned-spaceflight. Rendezvous 
technology has also evolved with small spacecraft development, such as the NASA’s 
Demonstration for Autonomous Rendezvous Technology (DART) [16], Air Eorce 
Research Eaboratory’s (AERE) Experimental Satellite Systems (XSS)-IO and XSS-11 
[17][18], Naval Research Eaboratory’s (NRL) Spacecraft for the Universal Modification 
of Orbits (SUMO) [19], and DARPA’s Orbital Express [20]. The SUMO program has 
evolved into the Front-end Robotics Enabling Near-term Demonstration (FREND) 
program, which maintains the objective of autonomous rendezvous and grapple 
operations. For a full discussion on the SUMO/FREND in context of the current status of 
autonomous rendezvous and capture missions, refer to Creamer [21]. Although, the most 
well known of these programs may be the XSS-11. The AERE Space Vehicle Directorate 
at Kirtland Air Force Base in New Mexico developed the XSS-11 in order to exhibit the 
ability for a small satellite to autonomously plan and rendezvous with passive and 
cooperative objects in EEO [22]. The use of micro-satellites to monitor, inspect, service, 
repair, and re-fuel larger spacecraft is a long term goal. The closest the XSS-11 
approached and maneuvered around another object in space was approximately 500 
meters. The XSS-11 used on-board laser range finders to measure the distance to target 
objects. Most of the XSS-11 flight is being conducted manually with autonomous 
planners running in the background. In addition, DARPA’s Orbital Express Advanced 
Technology Demonstration Program is intended to validate the technology and 
techniques for on-orbit refueling and reconfiguration of two satellites. The mission, 
which launched in late 2006, is intended to perform seven autonomous rendezvous and 
capture scenarios [23]. These will include component exchange and propellant transfer 
events. There is also research to apply rendezvous technology to smaller spacecraft, such 
as DARPA’s Tiny, Independent, Coordinating Spacecraft (TICS) program [24]. These 
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small spacecraft programs are leading the way for advanced autonomous close proximity 
operations by supporting the development of enabling technologies. 

There are several research groups contributing to the development of autonomous 
formation flight and docking control algorithms. Examples of recent progress in the field 
include such projects as Massachusetts Institute of Technology's (MIT) Space Systems 
Laboratory’s Synchronized Position Hold Engage Reorient Experimental Satellites 
(SPHERES) [25] [26] and AFRL TechSat 21 projects. The TechSat 21 program, which 
was cancelled in 2003, was intended to develop the concept of using clusters or 
formations of collaborative small satellites to fulfill complex missions, such as distributed 
aperture remote sensing, geo-location, of moving ground targets. As with most 
fractionated spacecraft architectures, the key motivating design concepts were that 
distributed mission architecture is more tolerant to damage and can also reduce the 
overall system cost. The SPHERES are intended to be used as a test bed for formation 
flight and reconfiguration, as well as autonomous rendezvous and docking technologies 
[26]. The first SPHERES satellite reached the ISS in May 2006 and has begun testing 
scenarios. These major research activities emphasize the challenge and ongoing 
requirement of developing a multiple spacecraft control algorithm. 

Maintaining tight formations for interferometry and other cooperative missions 
for long durations, is extremely taxing on spacecraft fuel. Also, thruster firings result in 
out gassing of propellant in close proximity may cause problems for spacecraft sensors 
and payloads. So, one way to ensure proper positioning and attitude relative to multiple 
spacecraft is to actually connect them in some manner. The connection of spacecraft on 
orbit is referred to as docking. Efficient and safe docking is the ultimate close proximity 
operational goal. The ability to dock spacecraft, in a safe, robust, and autonomous 
manner, is the cornerstone in being able to re-supply, service, upgrade, and reconfigure 
spacecraft while in orbit. 

1. Homogenous and Heterogeneous Spacecraft 

Any given group of spacecraft may have characteristics that make them similar or 
different. The likeness among vehicles is of significant consideration in the coordinated 
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control of a group of spacecraft. Spacecraft are termed homogenous, if the spacecraft are 
identical, or heterogeneous, if the spaeeeraft are different. The likeness evaluation of 
spaeecrafts may be based on a large variety of parameters, varying from physieal 
oharaeteristies (e.g., mass, volume, and strueture/eonfiguration) to the spaeeerafts 
capabilities (e.g., communication, data handling, power, payloads, and control systems). 
Limited types of launch vehicle capabilities and spacecraft bus designs give an upper 
bound to the mass and volume of orbital spacecraft. For instance, the static launch 
envelope inside the launch vehicle fairing limits the volume of a spacecraft [27]. 
Increased modularity of systems allows for a larger variation in payloads, sensors, and 
configurations. This research will be based on homogenous spacecraft, including similar 
control system sensors and aetuators. However, dissimilar size and shape spaeeeraft may 
be eonsidered as long as the dynamics of each spacecraft is properly modeled. Exploring 
the varying perturbations dynamics of multiple heterogeneous spacecraft is beyond the 
scope of this work. In this research, all commanded spacecraft are assumed to use the 
same basic control scheme and be equipped with sensors and actuators which offer the 
same level of precision. 

2. Orbiting Spacecraft Mission Phases 

The multiple spacecraft control algorithm used depends on the phase of spacecraft 
operations. In order to distinguish between phases of spaeeeraft operations based on 
physieal proximity, it is helpful to adapt some terminology from missile engineering. In 
missile intereeption, four fundamental stages of flight have been defined [28]. These 
stages are eommonly referred to as launeh, mideourse, terminal, and endgame stages of 
flight. They can be extended in spaeeeraft phase proximity operations, referred to as 
launch, midcourse, rendezvous, and docking. Here are the four spacecraft proximity 
phases: 

1. Launch phase ends after the satellite separates from the launch vehicle upper 
stage (booster) and it is in operational orbit. All large orbital maneuvers are 
performed during the launch phase. 
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2. Midcourse phase begins when the spaeeeraft has been stabilized into its 
operational orbit and it ean perform station keeping [29]. Small orbital 
eorreetions may be performed during this mideourse phase. 

3. Rendezvous, or terminal, phase is the when the spaeeeraft eonverges to a 
eommon point in its operational orbit. This rendezvous phase begins when 
the spaeeeraft reeeives its rendezvous eommand to translate toward a point in 
spaee. This spaee may be oeeupied by another spaeeeraft or be an empty 
spaee whieh the spaeeeraft will move into and oeeupy. If a group of 
spaeeeraft are eommanded to rendezvous to empty spaees relative to eaeh 
other they will form a spaeeeraft formation. On the other hand, if a group of 
spaeeeraft are eonverging on the same loeation in spaee they are eonverging 
for booking. 

4. Booking, or endgame, phase begins when spaeeeraft on-board proximity 
sensor aoquire the target looation/vehiole. Spaeeeraft booking is a very 
preoise maneuver, sinoe unoertainties in guidanoe and attitude need to be 
oorreoted quiokly and effeotively during this engagement soenario. For this 
researoh, short duration, olose proximity operations will be limited to 1.0 km 
of separation between spaeeeraft in nearly oiroular orbits. These assumptions 
are oonsistent with the staging required to get oooperative spaeeeraft into the 
same spatial region. Long term fuel optimal formation keeping and larger 
orbital eooentrioity variations are oonsidered separate issues whieh are not 
addressed in this researoh. 

3. Multiple Spacecraft Groups 

Before addressing the oontrol of multiple spaeeeraft, it is useful to define the 
differenoes between oonstellations, formations, and swarms/olusters. Constellations are 
groups of spaoeorafts in relative motion, or orbit(s), but their positions and attitudes are 
not dynamioally ooupled in any way [30]. This means that a ohange in position or 
velooity of one spaeeeraft does not impaot the others. For instanee, the GPS eonstellation 
is made up of at least 24 satellites with four satellites in six different GEO planes. Eaeh 
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satellite is maintained in its orbit independently via ground commanding. Maintaining 
the orbital position of the spacecraft, referred to as station keeping, is required due to 
several perturbation influences (e.g.. Sun and Moon gravitational fields, variations in 
Earth's gravitational field, aerodynamic drag). Constellations are usually maintained via 
centralized ground control. 

Formations are groups of spacecrafts which are dynamically coupled through a 
control law [30]. The motion of one spacecraft can give relative state information about 
at least one other spacecraft. Formation flight control of multiple small spacecraft is the 
task of maintaining the spatial formation, via control of the motion of the individual 
spacecraft in order to maintain the overall formation shape. The control law typically 
uses the position and velocity of one spacecraft to command another spacecraft [30]. The 
lead spacecraft moves along a commanded path and the following spacecraft(s) maintain 
a relative position, attitude, and velocity with respect to the leader. This is a common 
concept in leader/follower tracking schemes [31], where the relative motion is like a flock 
of geese in flight. Collisions are generally avoided through strict control of the following 
spacecraft's motion. For spacecraft formation control a virtual structure control approach 
is commonly used [11]. In spacecraft orbital terms, the lead spacecraft is referred to as 
the “Target” and the tracking spacecraft is called the “Chaser.” The target motion may be 
represented by an imaginary spacecraft. The idea of an imaginary target leads better 
understanding of a swarm/cluster. 

A swarm/cluster is a group of spacecrafts that move in concert with one another, 
but without strict control of relative positions, attitudes, or velocities. In the idea of a bee 
swarm, an outside observer can see the relatively smooth motion of an entire swarm but 
can not determine the specific relationship between any two spacecrafts. Think of the 
center of the swarm as an imaginary target and all of the spacecrafts as Chasers which 
only stay within a specific range. A tight swarm would be represented by a small 
acceptable distance between the Target and Chasers. However, the specific position, 
attitude, and velocity of each Chaser spacecraft are not centrally controlled. Each 
spacecraft autonomously manages its own motion and is only influence by its range from 
the target and the need not to collide with other spacecrafts. 
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The desire is to offer a control algorithm which bridges the gap between single 
cost optimal fuel trajectory tracking of rigid formation and the emergent behavior of 
swarms. Using APT based control, allows for incorporation of collision avoidance 
directly into the close proximity control algorithm. The goal position of each spacecraft 
is explicit and the general spacecraft path will be in the predictable direction of the goal. 
The obstacle potential functions can be defined and used to determine navigation paths 
which are robust enough to allow for both convergence and collision avoidance. For a 
detailed development of APF control refer to Chapter V. 
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III. SINGLE SPACECRAFT KINEMATICS AND DYNAMICS 

MODELING 


A. INTRODUCTION 

The first computational step of developing a control algorithm is to establish the 
system model. For this research, the fundamental system is a high fidelity six DOF 
spacecraft orbiting the Earth. This chapter introduces the dynamics and kinematics 
model used through the rest of the dissertation. The orbital perturbations included are 
fourth order harmonics in the Earth gravitational potential field, atmospheric drag on the 
spacecraft, third-body (Sun and Moon) forces, solar-radiation pressure, and mass 
variation due to thruster firings. Refer to Chapter X. for detailed description of the 
spacecraft model validation technique. 

A spacecraft orbit is usually expressed relative to a right-hand inertial (X, Y, Z) 
coordinate system with its origin at the center of the Earth and the center of the Earth in 
the plane of the spacecraft’s orbit. This reference frame is called Earth Centered Inertial 
(ECI), shown in Eigure 3.1. The X-axis points toward the vernal equinox, the Y-axis is 
90° counterclockwise from the X-axis in the equatorial plane, and the Z-axis extends 
through the North Pole [32]. Eor scaling reference, it is worth noting that the radius of 
the Earth spheroid is approximately, = 6,378.1 km. 
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Figure 3.1 Earth with ECI Coordinate System. 


A spaceeraft orbit is determined by the gravitational forees acting upon it. The 
spacecraft orbital dynamics, which is described in Chapter lll.B, can be simplified as the 
solution to a two body problem; refer to Chapter lll.B. 1. The description of a spacecraft's 
orbit can be characterized by a finite number of parameters, which are discussed in 
Chapter lll.B.2. However, the relative orbital position and velocity of a spacecraft 
depends on the reference frame of interest. Eor instance, in order to describe the 
spacecraft position with respect to the Earth, an inertial frame like that described above 
would be sufficient. However, there is also a need to describe one spacecraft position 
and velocity relative to another. Therefore, additional reference frames need to be 
defined; refer to Chapter lll.C. With objects in multi-dimensional space, there is not only 
a need to define a position and velocity in space, but also a need to define the objects 
orientation. The attitude dynamic of the spacecraft is discussed in Chapter lll.D. There 
are other forces and perturbations which influence the spacecraft orbit; these are 
discussed in Chapter lll.E. Einally, the overall characteristics of the six DOF spacecraft 
model used in this research are outlined in Chapter lll.F. 


18 















B. KEPLERIAN ORBITAL DYNAMICS 


I. Two Body Problem 

The fundamental description of a spacecraft orbit is determined by the 
gravitational force of the central body. In this document the Earth is assumed to be the 
primary body of interest and any spacecraft are in orbit around it. Simple two body 
orbital motion is described by Newton's Law of Universal Gravitation [33]. 



Gmm,( 




(3.1) 


where is the force of Earth's gravity acting upon the spacecraft, G is the Universal 
Gravitational Constant (6.673x10 5 ^^ fcg '), is the mass of the Earth 

(5.9733328x10^"^ kg), is the mass of the spacecraft, r is the relative distance vector 

from the center of mass (COM) of the Earth to spacecraft, and r is the Euclidean (2- 
norm) distance of the spacecraft from the COM of the Earth. A 3D position vector in 
Earth inertial reference frame, can be represented as 

r=r,X + r^Y + r^Z (3.2) 

The Euclidean (2-norm) of a 3D vector is 

r = 4{r,f+{r,f+{r,f=\r\ (3.3) 

The assumptions with this simple two body problem are that the Earth and spacecraft can 
be modeled as point masses, the gravitational field is symmetric, and no other external 
forces are acting on the Earth or spacecraft. Additional orbital forces and perturbations 
will be considered as the dynamic model is further developed in Chapter III.E. 

The left side of equation (3.1) can be expanded with Newton's Second Law [33], 

F = ma (3.4) 

where F is the force resulting from the mass (m ) and acceleration (a ) of an object. The 
relative acceleration of the spacecraft can be determined by subtracting the acceleration 
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of the Earth from the acceleration of the spacecraft, and substituting equation (3.1) and 
equation (3.4) results in 


a = - 








(3.5) 


The relative acceleration is the second derivative of the relative position vector j and 

the term on the far right is the normalized position vector. Next, make use of the 
gravitational parameter, defined as 

^ = Gm^ =398600.4418x10'(3.6) 


The value of // is based on modeling of the Earth, such as the World Geodetic System 
(WGS-84) [32]. Eor most applications, the mass of the spacecraft is negligible compared 
to the mass of the earth, so the two body motion equation simplifies as follows 





\rj 


= 0 


(3.7) 


This is the core equation for determining spacecraft position. Given the initial values of 
the relative position and velocity of the spacecraft, the orbit can be propagated using 
numerical analysis. The integration of Equation (3.7) leads to the classical Keplerian 
orbits, whose shape can vary according to the spacecraft initial position and velocity. In 
particular closed (circle, ellipse) and open (parabola, hyperbola) orbits are the possible 
trajectories for the spacecraft. 


2. Orbital Elements 

Any spacecraft specific orbit position can be also represented by six classical, or 
Keplerian, orbital elements, which can substitute the Cartesian coordinates. The six 
classical orbital elements are the following: 

1. Semi-major axis (a) determines the size of the orbital ellipse as 

r +r 

a = -—^ (3.8) 
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Orbits are typically referred to in categories that relate to general radial 
size of their orbit. Low Earth orbits (LEO) are typically around an altitude 
of 1,000 km. Middle Earth orbits (MEO) are typically around an altitude 
of 10,000 km. Geosynchronous orbits (GEO) occupy an altitude of 36,000 
km, and remain in a relatively stationary position relative to the Earth, 
depending on inclination. A special case of this type of orbit is the 
geostationary orbit, which remains in a stationary position over the Earth’s 
equator, due to an inclination of zero degrees. Using the semi-major axis, 
the period {Pj), or duration, of the spacecraft’s orbit can be determined as 



(3.9) 


2. Eccentricity (e) determines the shape of the orbit, such that when e = 0 
the orbit is circular. As the eccentricity increases from 0 up to 1, the shape 
becomes less circular and more ellipse-like. Eor e =1 the orbit is a 
parabola and for e>l the orbit is a hyperbola. 

3. Inclination (i) is the tilt angle relative to the Earth’s equator. If i = 0 
then the orbit is in plane with the Earth’s equator. 

4. Longitude of ascending node (Q) (also referred to as the right ascension 
of the ascending node) is the relative swivel angle of the orbital ellipse. 
This angle is measured from the inertial X-axis to the line defining the 
point where the orbit crosses the equator moving from the south to the 
north. 


5. Argument of perigee (w) is the angle between the ascending node and the 
closest radial distance (lowest altitude) from the orbit to the Earth 
(perigee). 

6. True anomaly (y) is the angle from the perigee to the actual spacecraft 
location, thus it varies throughout the orbit. This can be thought of as the 
phasing of the spacecraft in the orbital ellipse, or circle. 
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These six orbital parameters are useful for visualizing the spacecraft's orbit, as 
shown in 0. 


Z 



Figure 3.2 Classical Orbital Parameters. 

C. REFERENCE FRAMES AND TRANSFORMATIONS 

When discussing positions, velocities, and attitudes of orbiting spacecraft, one 
must have a suitable reference system defined. Although useful, the classical orbital 
elements are not the only way to compute orbital propagation. In some instances, the 
spacecraft position (f ) and velocity (v ) in the ECI frame may be used for translational 
computations. The use of the right-handed ECI frame is intuitively satisfying since it is 
centered upon the primary object of interest, the Earth. However, the ECI coordinate is 
independent of the spacecraft and its orbit and does not meet all of our needs. Eor 
instance, on-board spacecraft sensors are not centered in the ECI frame. 
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1. Earth and Body Centered Reference Frames 

There are several coordinate systems which are based on the plane of the 
spacecraft orbit. These typically use the spacecraft orbit as the fundamental plane, and 
are denoted by the first two axes in a right-handed coordinate system. As an example, 
one such system with uses characteristics of both the ECI and classical orbital elements is 
the perifocal coordinate system (P,Q,W), as shown in Figure 3.3. In the perifocal 
coordinate system, the Earth is the origin, the P-axis points toward the orbital perigee, the 
Q-axis is 90° ahead in the orbital direction, and the W-axis is normal to the plane of 
orbit. The P-axis and Q-axis are in the orbital plane and orientated based on the 
argument of perigee. If the eccentricity changes and the perigee rotates, then the PQW 
coordinate system will also change. The inclination tilts the PQW coordinate system by 
tilting the orbital plane. Although, useful in some situations, the PQW coordinate system 
is not always well defined. As for all systems based on classical orbital parameters, 
circular orbits in the equatorial plane require special rules. This is due to argument of 
perigee and semi-major axis not being defined for a circular orbit, and the longitude of 
ascending node not being defined for an equatorial orbit. Since the ECI and PQW 
coordinate systems centered at the same origin, the transformation from one system to the 
other is purely a matter of rotation. 


z 



Figure 3.3 Perifocal and ECI Coordinate Systems. 
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There are numerous situations when a referenee frame eentered on the spaeeeraft 
or its relative motion is needed. Eaeh spaeeeraft, assumed to be a rigid-body, has a body 
fixed-fixed referenee frame eentered at the spaeeeraft's eenter of mass. This spaeeeraft 
body frame (Xb, Yb, Zb) is typieally aligned with the prineiple axis of inertia. The 
spaeeeraft body frame is not inertial and is free to rotate as the spaeeeraft rotates along 
any of its three axes. The rotation around the body axes are represented by three angles: 
roll {(p), piteh {0), and yaw as shown in Figure 3.4. This frame is useful when 
determining the spaeeeraft attitude and rotation rates, similar to the referenee frame used 
for aireraft. 


Zg axis 



Figure 3.4 Body-fixed Coordinate System. 

Of partieular interest is the spaeeeraft eoordinate system (R,S,W) whieh is used to 
determine relative motion between objeets in orbit. The spaeeeraft eoordinate system is 
aligned with the R-axis along the radial direetion from the Earth to the spaeeeraft, the S- 
axis is along the direetion of the spaeeeraft translational traek, and the W-axis is eross- 
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track (normal) to the orbital plane. One possible orientation of the spacecraft, , as shown 
in Figure 3.5, is with the R-axis aligned opposite of body yaw, S-axis aligned with body 
roll, and W-axis aligned opposite the body piteh axis. In Figure 3.5, the spaeecraft 
coordinate system W-axis points directly up from the page and the body-fixed Yb axis 
points directly down into the page. As with the ECI and PQW coordinate systems, the 
body-fixed and RSW coordinate systems are centered at the same origin. Therefore, 
transformation from body-fixed to RSW is purely a matter of rotation. It is worth noting 
that the velocity vector is only aligned with the S-axis when the orbit is perfectly circular, 
and at the apogee and perigee of elliptieal orbits. The angular difference from the loeal 
horizontal (line perpendicular to the radial vector) and the velocity vector and is usually 
referred to as the flight path angle. 


Radial Line 
from Earth 



S and X| 


e 


W (toward reader) 

Yg (away from reader) 


R and -Z, 


B 


Orbital Path 


Figure 3.5 Body-fixed and RSW Coordinate Systems. 


The RSW and body-fixed referenee frames are not inertial. The spaeecraft 
coordinate system orbits along the orbital path in the ECI referenee frame, as shown in 
Eigure 3.6. 
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Figure 3.6 ECI and Spacecraft Body Coordinates. 

For spatial ranging and positioning information, each spacecraft will have on¬ 
board sensors and receivers. These sensors are likely to include Global Positioning 
System (GPS) receivers and directional range finders (e.g., sonar, laser, or infra-red) to 
determine distances away from the desired location and from other objects. Most payload 
sensors, with the exception of GPS receiver data, give information relative to their 
position on-board the spacecraft. Figure 3.7 shows an elliptical, equatorial orbit with the 
spacecraft coordinate W-axis and ECI Z-axis both pointing directly out of the page. In a 
2D or 3D space, ris the actual geometric distance between objects in the space, as 
computed in Equation (3.3). In order to transform data from the RSW to the ECI 
coordinate system both a rotation and translation is required. 
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Figure 3.7 Euclidean Distance. 


2. Transformations 

Transforming between two different reference frames with the same origin 
requires a rotation. For reference frames that do not have the same origin, the axis can be 
aligned by rotation and then translated with vector addition. For Cartesian 3D systems 
there are three degrees of freedom for angular displacement and one degree of freedom 
along each axis, resulting in nine element matrix with six relationships between them. 
The coordinate transformation matrix ( )> ^l^o referred to as a Direction Cosine 

Matrix (DCM) [34], is represented by 


c{0)c{y/) 


C 

'^DCM 


s{(p)s{e)c{ii/)-c{(p)s{y/) 
c{(p)s{9)c{y/) + s{(p)s{y/) 


c{0)s{y/) 

s{(p)s{9)s{y/) + c{(p)c{y/) 
c{(p)s{9)s{y/)-s{(p)c{y/) 


-s(9) 

s((p)c(9) 

c((p)c(9) 


(3.10) 


where the angles of yaw ( ^ ), pitch (9), and roll (^) correspond to the rotation described 
in the body-fixed reference frame, c(—) is the cosine function, and s(—) is the sine 
function. The order of angular rotation is important, since different order of rotation will 
result in a different transformation matrix. The , shown in equation (3.10), is for a 
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yaw, pitch, roll sequence. There are 11 other possible sequenees and resulting 
transformation matrixes. The matrix is useful for physieal insight of the 

transformation, but it is not effieient for eomputation due to the trigonometrie funetions. 

The most used parameterization for spaeeeraft attitude transformations are 
quaternions. The definition of quaternion is based on Euler's theorem whieh states that 
the general displaeement of a rigid body with one fixed point is a rotation about a fixed 
axis. This axis is ealled the eigenaxis or Euler axis The angle of rotation is ealled 

the Euler angle or the prineipal Euler angle (a). The Euler axis is the eigenveetor of the 
rotation matrix assoeiated with the eigenvalue 1. Thus, every rotation matrix has an 
eigenvalue that is equal to positive one [35]. This faet justifies the term eigenaxis for the 
Euler axis veetor, denoted by 


e ■ = 

axis 


(3.11) 


ft is useful to eompute Euler angles from a given rotation matrix, sinee there is a 
need to be able to eompute and a . The Euler parameter set, also known as a 
quaternion, is a four-parameter set with some eomputational effieieney and singularity 
advantages over the Euler angles set. Although useful, it may be diffieult to visualize the 
physieal meaning of the quaternion for most applieations. The quaternion veetor (q) is a 
3x1 [34], represented by 


q = 


gj sin(a / 2) 
sin(a / 2) 
^3 sin(a / 2) 


(3.12) 


The quaternion ((5 f)isa4 x 1 algebraie veetor [34], with a sealar eomponent as 
the fourth eomponent 

sin(a / 2) 

^2 sin(a / 2) 

q = 

^3 sin(a / 2) 

_eos(a/2) 
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(3.13) 



The spacecraft coordinate frame is not inertial, therefore it is moving with respect 
to the ECI frame. The kinematics of a non-inertial frame can be represented by the 
following equation for the velocity of the spacecraft with respect to the ECI coordinate 
system 

VEci=f + ^BN^r (3.14) 

The acceleration of the spacecraft with respect to the ECI coordinate system, can be 
determined by differentiating equation (3.14), 

^£C/ ^ ^ + ^ ^) + i^BN ^ ^) + ^BN ^ i^BN ^ ^) (3.15) 

These results can be generalized with respect to any two reference frames, and will be 
used to derive the equations for the relative motion between two spacecraft, refer to 
Chapter IV. 


D. ATTITUDE DYNAMICS 


I. Spacecraft Rotational Dynamics 

Applying Newton's Second Eaw to rotational dynamics, the torque on a system is 
defined by 


T=H=Ja 

The angular momentum is described by 

H = Jd> 


The Inertia matrix of a 3D system is a three by three matrix. 


J J 

XX xy 

J J 

yx yy 

■^zx ^ zy 


J 

J 

J 


xz 

yz 

zz 


(3.16) 


(3.17) 


(3.18) 


where the diagonal terms are Moments of Inertia and the off-diagonal terms are Products 
of Inertia. The inertia matrix is positive definite and symmetrical, with the Moments of 
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Inertia always positive or zero. The primary axes of the body can be defined along the 
principal axis, which are defined by the minimum and maximum moments of inertia. 

The Euler Equation for rotation follows the form of equation (3.14), and describes 
the torque in each axis [29] as 

ii + axH=f (3.19) 

Eor the system with the origin of the body frame at the center of mass, the Euler 
Equations [29] can be described by 

(a) 

(b) (3.20) 

(c) 

These three first-order differential equations are in components along the body 
axes of the angular velocity of the body frame with respect to the inertia frame. They are 
non-linear and coupled, so they are usually solved by numerical integration. The angular 
rate of the spacecraft body with respect to the initial frame is given by 

^BN ~ ^BO '^^ON (3-21) 

where is the angular velocity of the spacecraft body frame with respect to the orbital 
reference frame and is the angular velocity of the orbital reference frame with 

respect to the inertial reference frame. Eor a circular orbit with small angles and angular 
rates, the angular velocity of the spacecraft with respect to inertial is 


^ON ^Z 

(a) 


II 

1 

f 

(b) 

(3.22) 

COx CC^ + COqj^ CCx 

(c) 



2. Gravity Gradient 

Gravitational forces acting upon an orbiting spacecraft depend on the mass 
distribution of the body. Since most spacecraft are not perfectly symmetrical, the 
gravitational forces on a spacecraft are not uniform and result in a disturbance torque 
around the spacecraft's center of mass. The gravity gradient tends to align the 
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spacecraft’s minimum axis of inertia with the radial vector from the Earth. For example, 
a pencil orbiting in space will tend to align with the point toward the Earth. The 
Gravitational Gradient (GG) torque can be expressed along the principle axes of the body 
frame as 


T 

^GCx 

T 

^GCy 

T 

^GGz 


\r^ J 

y j 
\r^ J 


(Jz 

(Jx 


^ y ) G G 

■^z)GG 

■^x)GG 


{a) 

(b) 

(c) 


(3.23) 


where c^, Cj, and are the direction cosines of the radius vector with respect to the 

spacecraft's body frame. The GG torque is perpendicular to the local vertical. Therefore, 
the GG torque only has components in the roll and pitch axis. In fact, if one of the 
principle axes is aligned with the local vertical, then the GG torques is null. The GG 
torque can be used as a simple control method in order to stabilize a spacecraft attitude. 


E. ORBITAL PERTURBATIONS 

Non uniform Earth mass distribution and non spherical Earth shape, atmospheric 
drag on the spacecraft, third-body (Sun and Moon) forces, and solar-radiation pressure 
act as disturbances on a spacecraft Keplerian orbit [32]. The significance of these forces 
often depends on spacecraft size, position and altitude. In this research, interest is in 
relatively small spacecraft (less than 400 kg) at low orbital altitudes (less than 3,000 km). 
These conditions determine which perturbation forces have the largest influence on the 
spacecraft orbit. For relative low latitudes, the Earth's oblateness and atmospheric drag 
are more dominant. Also of importance is the time duration of interest. Some of these 
force are of great significance over the course of a satellites lifetime, but negligible for 
the duration of a single orbit. Simulations will include perturbations when practical 
based on the time duration and spacecraft characteristics and orbit. 
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1. Non-Symmetrical Earth 

The Earth is not symmetrical in its overall shape and mass distribution. This 
results in the gravitational field also not being symmetric. The Earth's oblateness, or 
equatorial bulge, can be described by the zonal harmonics coefficients in conventional 
Eegendre polynomials, refer to [32] for detailed discussion. The first three zonal 
coefficient terms in the Eegendre polynomial are J2, J3, and J4, where 

12 = 1.08262668355x10'^ is the second-order zonal harmonic. The third-order zonal 

-6 

harmonic is 13 = 2.53265648533x10 and the fourth-order zonal harmonic is 
J4 = 1.08262668355 x 10'\ J2 is the equatorial bulge term which has the most significant 
effect on spacecraft orbits. J2 effects are often classified as short period oscillations, 
while J4 results in long period variations. 

In this work, the EGM-96 coefficients and WGS-84 reference shape are used for 
calculations, refer to [32] for further details. Various, detailed Earth models are available 
through applications, such as STK [36]. Higher order polynomials will not be discussed 
in this research, although may be included in imported dynamic models. 


2. Atmospheric Drag 

The density of particles in the Earth's atmosphere is variable and changes due to 
solar interaction and magnetic field influences. Particles in the altitude of the spacecraft's 
orbit act upon the body of the spacecraft and slow it down. Atmospheric drag can impart 
both a translational disturbance force and rotational (attitude) disturbance torque. The 
basic equation of aerodynamic drag [32] is 


a, =- 

drag 2 


f A ^ 


6 - A 
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V ^rel y 


(3.24) 


where Cj^ is the spacecraft drag coefficient, is the spacecraft velocity vector with 
respect to the atmosphere, is the spacecraft’s cross-sectional area normal to its velocity 
vector, is the spacecraft’s mass, and p is the atmospheric density. The velocity 
vector relative to the Earth’s rotating atmosphere is 
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(3.25) 


where is the angular rotation veetor of the Earth. 

Comparable to the gravity models, numerous atmospherie models are available 
for referenee. In this work, an exponential atmospherie density model based on the U. S. 
Standard Atmosphere (1976) was used [32], This model, adopted by the U. S. 
Committee on Extension to the Standard Atmosphere (COESA), should be suffieient for 
EEO orbits with orbital altitudes between 100 km and 1,000 km [32], 


3. Third Body Effects (Sun and Moon) 

The gravitation attraetion from other bodies, such as the Sun and Moon, also 
affects spacecraft orbital Keplerian motion. The three body equation of motion for the 
relative acceleration of a spacecraft in the Earth's inertial frame [32] is 
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(3.26) 


where is the relative distance from Earth to Moon, is the relative distance from 
Earth to Sun, is the relative distance from spacecraft to Moon, is the relative 

distance from spacecraft to Sun, jU^ = 4902.798882x ^ is the Gravitational 

Parameter of Moon, and = 13271.2428xl0'^m^j'^^ is the Gravitational Parameter of 

Sun. If the terms in equation (3.26), especially with respect to the Sun, are too small for 
numerical precision then the acceleration can be approximated using Taylor series 
expansion [32]. 


4. Solar-Radiation Pressure 

The radiation being emitted by the Sun exerts a force on the spacecraft, somewhat 
like the wind on a sail. The magnitude and direction of the solar-radiation force is 
dependent on several factors, such as the position of the Sun relative to the spacecraft, the 
attitude and shape of the spacecraft, the intensity of the solar-radiation and the reflectivity 
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of the spacecraft [32], Solar pressure can impart both a translational disturbance force 
and rotational (attitude) disturbance torque. For this research, the spacecraft is treated as 
a “black body,” which absorbs all radiation. 

In general, the force of solar pressure per unit area on a spacecraft orbiting Earth 
is 5^ =4.51x10A small spherical spacecraft at lower altitudes will be 

influence less by this solar-radiation force. There may even be durations of the orbit 
where the spacecraft may be in eclipse and does not experience this effect. Although not 
expanded here, more precise equations for the force can be developed using the 
reflectivity of the spacecraft and the exposed area of the spacecraft. 


5. Thrust 

Thrusting, and any other control actuator output, from the spacecraft can be 
considered as a perturbation to the orbit. The thrust is the result of a controller's 
command to the propulsion system, in order to attain a desired orbit or attitude. The 
force of the thrust is 


F =I 

thrust sp 


^ dm 
V dt 


(3.27) 


where is the force of spacecraft motor, 7^^ is the specific impulse of the engine, 

dffi 

and - is the spacecraft motor mass flow rate. The change in velocity, Av, of the 

dt 

spacecraft due to the thruster firing is 


Lp In 


m„ 




(3.28) 


where g = 9.80665 m/s^ is the gravitational acceleration at the Earth's surface, and is 

the mass of the burnt propellant. The minimum and maximum Av are determined based 
on the duration of the thruster firing. Eor a given thruster force, the Av can be 
approximated as 
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F At 

= — ^2 29 ) 

Thruster sizing for a spacecraft depends on mission and orbit. Station keeping, 
drag compensation and three axis stabilization requirements all contribute to the total Av . 
In this work, each thruster is assumed to be hydrazine fueled with the compiled 
parameters listed in Table 3.1 [37]. The thruster(s) are assumed to be directional and 
mapped to deliver the commanded Av in the body frame of the spacecraft. 


Thruster Parameter 

Thruster Range 

Simulation Baseline 

Propellant Type 

Hydrazine 

Monopropellant Hydrazine 

Force 

0.5 -2N 

IN 

Specific Impulse () 

50 - 200 sec 

200 sec 

Steady State Impulse 

50 - 200 sec 

200 sec 

Min pulse duration 

0.01 - 0.5 sec 

0.05 sec 

Min Av 

0.01 - 5 mm/s 

0.0005 m/s 

Max total Av 

10 - 800 m/s 

60 m/s 

Mass 

3-7% of Spacecraft Mass 

3% of Spacecraft Mass 


Table 3.1 General Thruster Parameters. 


An on-off thrust profile is never ideal, due to variations in ignition, tail-off, mass 
flow rate, and specific impulse [32]. These variations can be managed by applying thrust 
modulation methods. The pulse width (duration of fire) and pulse frequency (rate of fire) 
can be used to limit the hysteresis, or dead zone lag. Using pulse width pulse frequency 
(PW-PF) modulation an approximately linear response can be achieved from on-off 
thruster. 


F. SPACECRAFT CONTROL SYSTEM 

The spacecraft control system, referred to as the Attitude and Orbital Control 
Subsystem (AOCS) [27], is based on the spacecraft’s kinematics and dynamics. 
Determining the actual position and velocity of the spacecraft in an orbit is the function 
of the Navigation, Guidance, and Control (NGC) subsystem [27]; whereas, determining 
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attitude is the function of the Attitude Determination Control Subsystem (ADCS) [27]. A 
typical control system block diagram is shown in Figure 3.8. If the attitude of the 
spacecraft is not important then the spacecraft can be treated as a point, or sphere, mass. 
However, most payloads require a specific pointing direction for their sensors, so the 
attitude of the spacecraft must be maintained. As the spacecraft location converges, the 
attitude of the spacecraft becomes more significant. For instance, if two spacecraft are to 
dock then their docking mechanisms must be aligned. The station keeping and attitude 
control (pointing) accuracy of spacecraft often vary significantly and are dependent on 
their sensors and actuators. Spacecraft attitude sensors include, earth. Sun, and star 
sensors, gyroscopes, magnetometers, and differential GPS receivers. Typical attitude 
actuators are divided into two categories; passive and active. Passive actuators include 
gravity gradient stabilization, spin stabilization, and dampers. Active actuators consist of 
thrusters, magnetic torquers, and momentum control devices (e.g., momentum wheels, 
reaction wheels, and control moment gyroscopes). 


Perturbations 

& 

Disturbances 



Figure 3.8 Control System (adapted from [27]). 

1. Translational Control 

Translational control is the major theme of this dissertation, with translational 
control achieved by position and velocity state feedback. The translational control is 
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discussed in detail in Chapters V, VI, and VII. Chapter V is an overview of spaeeeraft 
translational eontrol eoneepts. Chapter VI introduees the speeifie APF eontrol algorithm 
and Chapter VII diseusses the LQR/APF eontrol algorithm. 

Thruster firings, or any mass expelling system, ean be used to propel a spaeeeraft 
and guide its eenter of mass position and veloeity in spaee. The primary kinds of 
thrusters used on spaeeeraft are eold gas, hot gas and eleetrie. They eaeh have torque and 
speeifie impulse limitations, based on equation (3.28), as 

(3.30) 

mg 

The propulsion is primarily dependent on the mass flow rate. The advantage of 
using thrusters is that they ean be used to torque about any axis in any orbit. However, 
thrusters generally have on-off limitations whieh may result in sensitivity to ehattering, 
dead zones, and hysteresis issues. These issues are of partieular eoneem in the final 
stages of doeking, due to the need for preeise eontrol. These ehallenges ean be overeome 
by modulating the thruster firing duration and frequeney. But these modulation sehemes 
add to the overall eomplexity of the eontrol law being used. 

Also, various sizes and numbers of thrusters ean be positioned at desired loeations 
on the spaeeeraft. Smaller thrusters are used for fine adjustments and larger thrusters are 
used for larger maneuvers. In spaeeeraft motion thrusters are typieally used to impart a 
Av in the desired direetion and then another Av in the opposite direetion onee the 
spaeeeraft is elose to the desired spatial region. This seeond thruster firing serves as a 
braking maneuver, sueh as a ear approaehing a stop sign. For short duration firings, 
thruster firings ean be treated as generating instantaneous veloeity variation of the 
spaeeeraft, not ehanging its position. 

2. Attitude Control 

Attitude eontrol was aehieved via quaternion feedbaek eommands to three axis 
aligned reaetion wheels with magnetotorqers for momentum damping. The non-linear 
quaternion feedbaek eontrol law utilized is based upon on attitude error deseribed by the 
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quaternion and the angular rate (body with respeet to inertial). The eommanded torque 
due the quaternion error and angular rate is 


^CMD ~ ^BN (3.31) 

where and are the proportional and derivative feedbaek eontrol gains, respeetively. 

These gains were tuned in order to get a damped attitude response with an aetuation 
limited maximum torque. Of eourse, the maximum torque performanee is dependent on 
the attitude eontrol aetuator. Additionally, quaternion rotation logie is neeessary to 
ensure minimum rotations are eommanded. Sinee the quaternion represents a veetor in a 
4-D unit sphere with an axis and an angle, it is possible for two different quaternions to 
represent the same orientation. In partieular if the opposite eigenaxis is eonsidered the 
angle rotation is 360 degrees less, represented as 

r 

or q = 

eos(«/2)J [eos(« / 2) - 2;r 

The quaternion with the positive fourth element represents the shortest rotation, and its 
eomplement is the longest. Seleetion of quaternion with minimum angular rotation is 
generally desired for effieient attitude eontrol. 

Momentum Exehange Deviees (MED), sueh as, momentum wheels, reaetion 
wheels, and eontrol moment gyros, are eommonly used for providing eontrol torque. The 
general attitude dynamies for a spaeeeraft with a MED is 

j ^BN ^BN ^ ^ ^BN —T + (3.33) 

whieh follows from Equation (3.17) and Equation (3.19). The spaeeeraft torque due to 
the MED () is deseribed by 

TMED^-{^^BN^^h) + iyZh^ + [zh) + [j (3.34) 

where h is the angular momentum of the MED and Z represents the axis about whieh the 
MED spins. The mass of the spinning wheel aets as inertia to the spaeeeraft system and 
provides disturbanee rejeetion in the perpendieular axis. 



(3.32) 
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Each of the three primary types of MED work in a slightly different manner. The 
simple momentum-wheel is a fixed-axis wheel spinning at a large angular rate. This is 
the most simple of the MED, due to the eonstant mass, fixed axis, and reasonably 
eonstant spin rate. Eor instanee, a torque due to momentum wheel rotating about the 
spaeeeraft’s piteh axis is 
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(3.35) 


where the last two terms on the right side of Equation (3.34) are now zero. This ean be 
simplified for a single momentum wheel, by substituting equation (3.22) into equation 
(3.35), as 


-{Kw 

h 

mw 

(h d )-(h CO a ) 

V mw X) \ mw mw z ) 


(3.36) 


The reaetion wheel is more eommonly used. The reaetion wheel is a fixed wheel 
that varies its spin to deliver desired torque. Three reaetion wheels are typieally aligned 
along the prineiple axes along with three magnetotorquers. The orientation of the 
reaetion wheel, Z , is a eonstant, so Z = 0 and the Z*h term in (3.34) ean be dropped. 
Eor a typieal spaeeeraft eonfiguration, with three reaetion wheel along eaeh axis and 
rotating in the same direetion; 


Z 




1 0 0 
0 1 0 
0 0 1 


(3.37) 


with the spin direetion defined to be with respeet to eaeh axis. Mass of the spaeeeraft is 
still eonsidered relatively eonstant, so / = 0 and the J * term ean be dropped. The 

angular rate of the reaetion wheel is variable and therefore h^O (different than 
momentum wheel). Therefore, the general eontrol law for the reaetion wheel ean be 
redueed to 
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which can be slightly expanded as follows 
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(3.39) 


It is useful to note that that the rate-of-ehange of the wheel momentum is as follows 

^RW ~ (~'^CMD ~ ^BN ^ ^RW * ^RW ) -^O) 

Equation (3.40) is used to deseribe the dynamies of the reaetion wheel. The eommanded 
torque is designed to saturate at the maximum torque available to the reaetion wheel 
along eaeh axis. 

The magnetotorquers may be used to de-saturate the reaetion wheels, as 
neeessary. The reaetion wheel momentum is multiplied by a gain, , in order to 

determine the torque required, , by the magnetotorquers eontrol law. 


^REQ 


^mt^w 


(3.41) 


This is used to determine the magnetie dipole momentum, M , of eaeh 

magnetotorquer, as 


M = 


R xT 

^ REQ 


(3.42) 


where is the magnetie field in the body frame. The maximum dipole moment of eaeh 

torquerod, M , is limited. Also, the piteh axis is not torquable at 11.7 degree orbital 
inelination sinee the spaeeeraft is at the magnetie equator. At the magnetie equator the 
piteh axis is aligned with the magnetie field and ean not ereate torque along that axis. 
With the limited dipole moment and the magnetie filed with respeet to the body frame, 
the de-saturation torque generated by the magnetotorquers is given as 

TDAMP=MxBg (3.43) 
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This torque tends to align the dipole with the magnetic field. 


Each spacecraft model will be assumed to be equipped with an attitude control 
system. The spacecraft attitude control actuator and damper reactions will be based on 
reasonable attitude control systems for small spacecraft. Standard spacecraft attitude and 
rate sensors will also be assumed, with typical statistical sensor noises. Development of 
specific sensors and actuators is beyond the scope of this research. 

It is worth noting that although thrusters are the most common actuators to control 
spacecraft translational motion they can also be used for external spacecraft torque. The 
thrusters are usually positioned at equal and opposite radial distances from the center of 
mass in order to create torque. Various sizes and number of thrusters can be used, but 
they usually range from four to twelve per spacecraft [29]. Thrusters can also be used to 
supplement momentum exchange devices, by desaturating spin rates and changing 
momentum by providing external spacecraft torque. 


G. HIGH FIDELITY SIX DOF SPACECRAFT MODEL 

The characteristics of each of the multiple spacecraft in the group are assumed to 
be the same unless explicitly stated in the simulation results. The spacecraft altitudes are 
limited to LEO orbits in the range of 300-2,000 km. The distance between the Target and 
Chaser spacecraft initial positions is limited to less than 1.0 km in RSW coordinates. The 
basic characteristics of the spacecraft considered in this research are listed in Table 3.2. 
These characteristics were selected to represent realistic and current operational 
capabilities of a small spacecraft. Rules for spacecraft and subsystem sizing were 
determined based on design rules of thumb [37]. Eor instance, the spacecraft volume was 
selected to be approximately 1/100 of the spacecraft mass [37]. The center of mass of 
the spacecraft is assumed to be located at the geometric center. Position, attitude, and 
ranging sensors are assumed to provide near ideal information. Translational motion is 
conducted via six cold gas pulsed thrusters with a maximum thrust of 1.0 V . The attitude 
actuators are three reaction wheels along each spacecraft axis, along with 
magnetotorquers to dump the reaction wheel momentum. Disturbance perturbations are 
adjusted as necessary for model and simulations development. The attitude actuator 
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specifications are loosely based on the Honeywell Model HR 0610 Reaction Wheel [38]. 
The translations thrust actuators were generalized from cold gas thruster performance and 
operating characteristics; refer to [37]. 


Spacecraft 

Physical 

Characteristics 

Length and Width 

1.0 m 

Height 

1.0 m 

Mass 

100 kg 

Moment of Inertia about X 

16.67 - 50 kg rev 

Moment of Inertia about Y 

16.67 - 50 kg m^ 

Moment of Inertia about Z 

16.67 - 50 kg rev 

Spacecraft 

Actuators 

Number of Thrusters 

1-6 

Mass of Propellant (% Mass) 

3-6% 

Max Thrust (each thrusters) 

LON 

Number of Reaction Wheels (RW) 

3 

RW Max Torque (each axis) 

0.055 Nm 

RW Max Angular Momentum 

4-12 Nms 

Initial Angular Rate of RW 

ORPM 

Inertia about RW spinning axis 

0.14258 Kg m^ 

Number of Magnetotorquers 

3 

Max dipole moment (each torquerod) 

100 A m^ 


Table 3.2 Main Characteristics of Chaser Spacecraft Simulators. 


Initially these spacecraft characteristics were used to model spacecraft using The 
MathWorks Incorporated products MATLAB and Simulink [39]. The spacecraft 
characteristics may be further expanded and refined as hardware-in-the loop structures 
are developed at the NFS Spacecraft Servicing and Robotics Laboratory. Mass and 
dimensions greatly affect both the perturbations on the spacecraft and the commanded 
actuation by the control algorithm. Variations in the moment of inertia for three basic 
spacecraft shapes, spherical, cylindrical, and cubic, are considered. Mass changes due to 
thruster firings are modeled with appropriate moment of inertia changes. In the absence 
of actual sensor readings and actuators responses, equivalently accurate input and outputs 
are provided into the simulations. For instance, attitude determination sensor reading, 
such as star trackers, can be simulated based on orbital propagation calculations. On the 
other hand, actual GPS reading in the laboratory may be used to determine spacecraft 
positioning on the Autonomous Docking and Spacecraft Servicing testbed. 
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IV. SPACECRAFT ORBITAL RELATIVE MOTION 


A. EQUATIONS OF RELATIVE MOTION 

As two spacecraft obit around the Earth their relative positions change based on 
the characteristics of their orbits. Two different orbits are shown in Figure 4.1. This 
view shows a circular and an elliptical equatorial orbit from the perspective of high above 
the North Pole. The difference in inertial position and velocity between spacecraft orbits, 
appears as dynamic relative motion between orbital spacecraft. 



E 


22 

X 

05 

I 

>- 



X 10 
1 . 




-0.5 


X-axis (meters) 


X 10 


Figure 4.1 Two Spacecraft Orbits in 2D. 


In order to establish the equations of motion between spacecraft consider one of 
the spacecraft as the primary and all of others as secondary. The primary spacecraft will 
be denoted as the Target (t). The secondary spacecraft will each be denoted as Chasers 
(c) with numerical designations to distinguish them apart. The Target and Chase 
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spacecraft reference is used during the rendezvous maneuver, when the Target spacecraft 
maintains its orbit and the Chase spacecraft maneuvers to the Target's position. 

The relative motion can be developed from the position vectors of the two 
spacecraft. From Figure 4.2, it is apparent that vector addition gives 

r^=r,+f (4.1) 

where i; and i; are the position vectors of the Target and Chase spacecraft in the ECI 
coordinates and f is the relative position of the Chase with respect to the Target. 


z 



Figure 4.2 Relative Distance. 


Using equation (3.14), with the Target spacecraft as the origin of the rotating 
reference frame, the relative translational acceleration of the Chase spacecraft is given by 

f^=f,+ 2(d)^^ X f) + (®^^ X f) + X x f) (4.2) 

In order to solve for the relative motion, the position vectors can be expressed in the 
Target spacecraft coordinates. The relative position vector is 

r=xR + yS + zW (4.3) 
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with the (x, y, z) adopted as standard Cartesian eoordinate notation. The position of the 
Target spaeeeraft is displaeed from the ECI eoordinate along its radial axis 

= —r^R (4.4) 

Substituting equations (4.3) and (4.4) into equation (4.1), the position veetor of the Chase 
spaeeeraft beeomes 

r^={x-r^)R + yS + zW (4.5) 

Considering the foree of gravity aeting on the spaeeeraft from the two-body problem 
from equation (3.7), the aeeeleration of the Target spaeeeraft is denoted as 



The aeeeleration of the Chase spaeeeraft is denoted as 

C = ~Sc + ^thrmt 


(4.6) 


(4.7) 


where is the aeeeleration due to eontrol aetions and is the aeeeleration of the 

Chase spaeeeraft due to the foree of gravity. Aeeeleration on the Chase spaeeeraft due to 
gravity ean be resolved into the Target spaeeeraft eoordinates [40], 


Sc 





(4.8) 


where is the veetor norm, as denoted in equation (3.3). Next, these resolved 

eomponents ean be substituted into equation (4.2) with the assumption that the angular 
veloeity of the of the Target spaeeeraft eoordinate system is 

= oM (4.9) 

This angular veloeity represents that the orbital rotation only oeeurs around the 
spaeeeraffs eoordinate axis perpendieular to the fundamental orbital plane. By 
performing the eross produets and simplifying, the aeeeleration equations of the Chase 
spaeeeraft along eaeh of the Target spaeeeraffs axis eomponents are 
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(b) (4.10) 

(c) 


This set of equations represents the nonlinear equations of relative motion [40]. 
Using the same orbits as shown in Figure 4.1, the relative motion of the Chase spacecraft 
elliptical orbit with respect to the Target spacecraft circular orbit is shown in Figure 4.3. 
The y-axis is the radial axis from the Earth and the x-axis is transverse position of the 
Chase spacecraft. The plot is orientated for a counterclockwise orbital rotation with the 
velocity vector pointing to the left. Notice that the motion of the Chase spacecraft drifts 
away from the initial (0, 0) position. Additional translational acceleration forces can be 
incorporated into the equations in the same manner as the gravitational and control 
forces; refer to Chapter Ill.E. 
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Several algorithms have been developed to implement these equations. They vary 
in aeeuraey based on the gravitational and referenee orbit assumptions. The most 
eommon is based on a spherical Earth and a circular reference orbit [41]. Although for 
long term formation flight, geometric methods which consider reference orbit eccentricity 
and differential gravitational perturbations may be more accurate [42]. However, as the 
relative motion methods increases in accuracy they also increase in complexity. Detailed 
analysis of the orbital eccentricity and perturbation effects can be evaluated in order to 
determine which relative motion theory should be applied to a relative motion of two 
spacecraft [43]. Although for this research pertaining to short duration proximity 
operations the Clohessy-Wiltshire equations are sufficient. 

B. THE CLOHESSY-WILTSHIRE EQUATIONS OF MOTION 

The nonlinear equations for relative motion, shown in Equation (4.10), may be 
simplified and linearized in order to provide closed form analytical solutions. The 
Clohessy-Wiltshire equations are the well known linearized relative motion equations for 
near circular orbits. Eirst, assume that the relative distance between the two spacecraft 
(r ) is much smaller than the orbital radius of the Target spacecraft (r,) 

Tj » r (4.11) 

This research assumes that the relative distance is within a few kilometers and that the 
orbital altitude is at least a few hundred kilometers. This leads to the following, 
linearized estimations [32] [40]; 
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(4.12) 


Substituting the equation (4.12) estimates into equation (4.10) yields the following 




x = i-g 


1-^ 

V r y 


+ a^ + g + 2(a)y) + d>y) + (co^x) > R 


y = - 


Z = i- 


V rj 
g- 

V r J 


+ a - 2((yi) - {d>x) + (co^y) [• S 


+ a. 


(a) 

(b) 

(c) 


(4.13) 


Next, assume that the Target spaceeraft is in a near circular orbit. This 
assumption yields a simple estimate of the orbital angular velocity, in radians per second, 
as 


(4.14) 

Based on this estimate the angular velocity is constant, therefore the angular acceleration 
is zero (<» = 0 ). These assumptions yield the simple linearized equations known as the 
on Clohessy-Wiltshire equations 


X - 2(coy) - 3(0)^ x) = 

(a) 


y + 2(cox) = a^. 

(b) 

(4.15) 

z + {(o^z) = 

(c) 



where the acceleration due to control actions is zero for freely orbiting spacecraft. 
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These equations ean be written in general state spaee form, x = Ax + Bu and 
y = Cx, with the following six states 

0 0 0 1 0 Olfxl To 0 0 

000 0 lOj 000 

000 0 Olz 000 

3co^ 0 0 0 2® Ox 100 

0 0 0 -2® 0 0 j 0 10 

0 0-®^ 0 0 0 [zj [0 0 1 

1 0 0 0 0 OlTx' 

0 1 0 0 0 0 J 

0 0 1 0 0 0 z 

(4.17) 

0 0 0 1 0 0 i 

0 0 0 0 1 0 j 

0 0 0 0 0 lj[z 

The nonlinear eoupling weakens as the range between spaeeeraft deereases and 
remains mueh smaller than the Target’s orbital radius [41], These state spaee equations 
ean be evaluated for eontrollability, observability, and stability. The linear system is 
eompletely eontrollable, sinee the eontrollability matrix, , is full rank. 

Q,=[b AB A5"-‘] 

and eompletely observable, sinee the observability matrix, , is also full rank. 

“C 
CA 

R = . 

o \ 

CA""‘ 

The unforeed relative dynamie system is unstable, sinee the six eigenvalues of the 
A matrix eonsists of two repeated roots at zero and two repeated imaginary eomplex 
pairs. For the zero relative veloeity initial eondition, the stability of the system ean be 
eoneeptually thought of as a two axis of stability along the eigenveetors of the zeros and 
planes deseribed by the eomplex eigenvalues related veetors. Relatively elose initial 
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(4.18) 


(4.19) 






Chase spacecraft positions along the S-axis, will stay within a bound and may be 
considered equilibrium points. Also, there is an axis slightly off of the S-axis with 
velocity changes along the same axis which may be considered an axis of equilibrium. 
However, the Chase spacecraft must meet the condition of being relatively close, within a 
few kilometers, and near the same circular orbit as the target spacecraft. The other 
possible equilibrium locations may be considered dynamically unstable due to their strict 
dependence on both position and velocity relationships. These purely imaginary 
eigenvalues explain the oscillatory behavior of spacecraft relative dynamics. The system 
stability is further addressed in Chapter VII.D. 

Using the same orbits as shown above, the linear relative motion of the Chase 
spacecraft’s elliptical orbit with respect to the Target spacecraft’s circular orbit is shown 
in Figure 4.4. Once again, the y-axis is the radial axis from the Earth and the x-axis is the 
transverse position of the Chase spacecraft, with the velocity vector to the left. Notice 
that the motion of the Chase spacecraft is now periodic and does not drift due to 
eccentricity of the Chase orbit. 
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Figure 4.4 Linear Relative Motion. 

Next, Figure 4.5 shows both the linear and nonlinear motions on the same plot. 
The Clohessy-Wiltshire equations traek the true relative motion for short time periods. 
The absolute relative differences between the nonlinear and linear equations for the orbits 
are shown in Figure 4.6. In this instance, the Clohessy-Wiltshire equations vary less than 
70 meters from true for the first 30 minutes of the orbital propagation. As a practical rule 
of thumb, any step duration longer than a quarter of an orbital period should be avoided 
due to the 90° rotation of the RSW coordinate system. The set of ordinary differential 
equations, in equation (4.15), have been used extensively in rendezvous and docking 
algorithms. They have been useful in analyzing trends for initial positions and velocities 
[32]. The relative motion equations of the spacecraft can be further simplified into a state 
transition matrix form; refer to Chapter IV. C. 
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Figure 4.6 














































C. STATE TRANSITION MATRIX FOR RELATIVE MOTION 


The unforced ordinary differential equations, in equation (4.15), can be solved 
using Laplace operators. The solutions of these Clohessy-Wiltshire equations are known 
as Hill's equations. The position solutions, based on constant initial values, [32] are 


^ X ^ 
'An 


x{t) - 1^—J sin(ry t) - + — J cos(ry t) + 


4xo + 


CO J 




6xq + 


4jo 



f 2i:„) 

sin(u;f) + 
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KCO j 


f 2x'' 

c(Mo}t)-[6coxQ+2>%)t+ y^H—2. 

V y 


z(t) = -Zq co^micot) H—^ cos(ry 0 

CO 


(a) 

(b) (4.20) 

(c) 


where the x^, y^, and Zq are the initial position values and x^, y^, and Zq are the initial 

velocity values for the Chase spacecraft with respect to the Target spacecraft's RSW 
coordinate system. The velocity solutions [32] are 


x{t) = Xq cos(® t) + (3(y Xq + 2 jg ) sm{cot) 

y{t) = [6coxQ +4yQ)cos((yt)-2ioSin((yt)-(6(yxQ +3yo) 

(a) 

(b) (4.21) 

zit) = -Zq cosm(cot) +—cos(cot) 

(c) 


CO 


The R and S coordinates are still closely coupled, and the W term is a simple harmonic 
oscillator. 

The position and velocity equations can be placed into general state-space form 

X(0 = OXo (4.22) 


where X(t) is the current state vector, O is the state transition matrix, and is the 

initial valued state vector. The state vector is composed of the position and velocity of 
the Chase spacecraft, defined as 


X =[x,y,z,x,y,zf =[f,ff 


The state transition matrix relating to this state vector is 


o = 


A 

C 


B 

D 


(4.23) 


(4.24) 
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where A, B, C, and D [40] are given as 


A = 


4-3eos((»0 0 0 

6(sin((»?)-o?) 1 0 

0 0 cos{cot) 


(a) 


B = 


1 2 

—sin((»?) — (l-eos((»0) 

CO CO 


0 


0 


0 


2 4 

— (eos((»0“l) — sin((»?)-3? 0 

00 00 


— sin(o?) 
00 


(b) 


C = 


3(»sin(o?) 0 0 

6co{cos{cot)-Y) 0 0 

0 0 -o(sin(o?) 


(c) 


D = 


cos{cot) 2sin(o?) 0 

-2sm{cot) -3 + 4cos{cot) 0 
0 0 cos{cot) 


(d) 


(4.25) 


This state transition matrix ean be used to effieiently ealeulate the initial and final 
veloeity of the Chase spaeecraft required to aehieve a desired position. Any required 
magnitude ehange in the Chase spaeeeraft's veloeity eorresponds to thruster firings, 
denoted as Av . All neeessary Av's ean be summed and used as a measure of the eost of 
the spaeeeraft's maneuver. The minimum total Av maneuver ean be determined as a 
baseline for evaluating eontroller eommand performance. This will be discussed in more 
detail in Chapter IV.D. 


D. OPTIMAL TWO IMPULSE RENDEZVOUS 

The co-planar two thrust maneuver, based on the Hohmann transfer concept, is 
the classic way to plan the rendezvous of two spacecraft. The Hohmann maneuver is 
mathematically accepted as the most fuel efficient transfer between two circular, and 
elliptical, planar orbits [44]. The maneuver, first suggested by Walter Hohmann in 1925, 
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involves two tangential thruster firings for start and stop of the transfer [45]. The two 
thruster firings are eaeh approximated as impulse ehanges in veloeity. This maneuver 
minimizes the Av required, but takes a longer time to eomplete the maneuver. The 
maneuver time is one half of the transfer orbit's period. There are eases where Hohmann 
transfers ean be used in series, ealled bi-elliptie transfers, in order to aehieve less Av 
[32]. The two impulse rendezvous maneuver can be used for fuel efficiency comparison, 
and the bi-elliptical maneuver will not be explored in this research. This orbital transfer 
maneuver is a classic example of a known fuel efficient result, which can be used to 
validate the high-fidelity six DOF spacecraft model. 

The two impulse maneuver is not physically possible or practical for spacecraft 
close proximity operation. The impulse velocity change is physically impossible for a 
spacecraft with non zero mass and momentum. The force required to gain the desired 
impulsive is not usually possible with the thrusters onboard spacecraft. Additionally, the 
second impulse occurs at the goal position and does not leave much safety margin. 
Position uncertainties and rapid convergence at the rendezvous of two spacecraft make 
the second impulse undesirable. Despite these limitations, the two impulse maneuver 
serves as a baseline for minimum fuel efficiency and maximum duration of a spacecraft 
maneuver between two points. Therefore, the two-impulse maneuver is useful, as the 
baseline, for comparing the performance of spacecraft maneuver control algorithms. 

The Hohmann maneuver thruster firings are predetermined to take place at a time 
corresponding to apogee of the initial and transfer orbit. Small change in the time of the 
thruster firing can result in drastic variation in the Av for a general point-to-point 
transfer. For evaluating the general point-to-point orbital transfer, rewrite the state- 
transition matrix for a final position (iy ) and velocity (iy ) at some final time 

{f^f =A{r,f+B{7j (a) 

(4.26) 

{-fff =C{rJ +D{rJ (b) 

where Fq = [-^o’3^o’^o] Fq = [-^o ’ Fo ’ the initial relative position and velocity. 

The starting relative velocity of the Target and Chase spacecraft is selected to be zero. 
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such that the Chase and Target veloeities are equal. To ealeulate the initial veloeity to 
intereept the Target position, equation (4.26) ean be rearranged 

{%f -A{fJ) (4.27) 

where the inverse of B is 

Aa>s,{n{co\)+l>o/\ -2(y(cos((yt)-l) ^ 

(-8+3 (y t sin(w t)+8cos(w t)) (-8+3 co t sin(w t)+8cos(w t)) 

26;(cos(6;t)-l) -6;sin(6;t) ^ (4.28) 

(-8+3 CO t sin(w t)+8cos(w t)) (-8+3 co t sin(w t)+8cos(w t)) 

0 0 —-— 

sin(6;t) 

For a given initial and desired final position, the final veloeity ean be ealeulated 
by substituting equation (4.27) into equation (4.26) to yield 

{yf =C(Fo)"+£>[5-‘((c)"-A(F„)")] (4.29) 

For the situation where the Chase spaeeeraft is starting and finishing at rest 
relative to the Target spaeeeraft, the total Av is determined by the initial and final 
veloeities. The initial ehange in veloeity is 

Avo=|Fo| (4.30) 

and the final ehange in veloeity is 

Av^ = -Vj (4.31) 

The final veloeity, from equation (4.29), is the aetual veloeity of the spaeeeraft at the 
final position. The final ehange in veloeity needs to be eounter aet this veloeity in order 
to stop the relative motion at the goal. The total ehange in veloeity of the maneuver is 
simply the sum of all of the two impulse ehanges in veloeity 

Av = Avg + Av^ (4.32) 

However, the value of Av depends on the time seleeted to perform the impulse 
thruster firings. Therefore eare must be taken in establishing the optimized Av , based on 
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a functional minimization algorithm which contains time dependent minima. This is 
better illustrated in the following two examples, whieh ean serve as a basis for 
eomparison. 

In the first ease, the Chase spaeeeraft is direetly ahead of the goal position. The 
Chase spaeeeraft is initially at [0, 70, 0] meters in the Target’s RSW eoordinate system 
and maneuvers to the Target’s position of [0, 0, 0] meters. For a data set sampled every 
seeond for approximately 20,000 seeonds, using the state transition equations, determined 
in Chapter IV.D., yields maneuver duration of 

C =16,2785 (4.33) 

and a total optimized Av of 

Av = 2.8636xl0^'m/5 (4.34) 

These results are eomparable to those in [40]. The orbital maneuver eould have been 
done more effieiently over a longer duration. 

In the seeond ease, the Chase spaeeeraft is offset from a goal position. The Chase 
spaeeeraft is initially at [50, -100, -50], maneuvers to the Target’s position of [0, 0, 0]. 
For a data set sampled every seeond for approximately 20,000 seeonds, the optimal time 
of transfer was determined to be 

C= 3,972 s (4.35) 

whieh yields a total optimized Av of 

Av = 0.19675 m/5 (4.36) 

These results are on the same order of magnitude as those in obtained by [40]. 

The optimal trajeetory path will vary if the final time is allowed to vary. Therefore, the 
optimal trajeetory deseribed here is a funetion of both the time and Av of the maneuver. 
Using different time eonstraints will result in different paths to the final position. This 
idea ean be applied to the situation where an obstaele is in the optimal path and a 
seeondary path needs to be seleeted in order to avoid eollision. This will result in an 
effieient, but not Av optimal, trajeetory. 
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Minimizing the Av based on the Clohessy-Wiltshire equations, will result in 
optimal trajectory for point-to-point transfer maneuver. However, this maneuver is not 
flexible and is compromised if any course correction is necessary [40]. Small changes in 
desired waypoints or timing of the impulse thrust can result in large trajectory variations 
and offsets from the desired position. However, these optimal solutions for minimum Av 
maneuvers can serve as useful baselines against which an artificial potential field based 
control algorithm can be evaluated. 
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V. CLOSE PROXIMITY SPACECRAFT CONTROL 


A. CONTROL ALGORITHM DISCUSSION 

There are several different eontrol strategies and sehemes that may yield 
satisfaetory results for a given seenario. They all have advantages and disadvantages. 
These strategies range from striet single eonstraint optimal eontrol to flexible adaptive 
eontrol. The optimal fuel eontrollers are useful for large spaeeeraft maneuvers whieh 
require minimum use of fuel. While, the flexible adaptive eontrollers are applieable for 
spaeeeraft in deep spaee that are likely to experienee large eommunications lags, sensor 
failures, and unknown situations. By evaluating the disadvantages of eaeh, it ean be 
determined that neither would be ideal for close proximity spacecraft operations. 
Optimal path planners have the disadvantage of being too slow and computational 
intensive for real-time application in the presence of obstacles [19]. While, adaptive 
schemes tend to exhibit more arbitrary motion and require hierarchical structures to 
maintain fuel efficiency. Although, previous consideration of artificial potential fields 
for on-orbit assembly [40], motivates the development and application of a potential 
function based controller for multiple spacecraft in close proximity operations. 

Optimal control is typically based on minimizing a single cost function, which 
may consist of many diverse parameters. Various parameters, including distance, time, 
energy and path/sensor mapping can be optimized for path planning. However, only one 
of these parameters can truly be optimal for any particular cost function. All others 
parameters are treated as constraints. As multiple parameters are included in a cost 
function the range of possible solutions exhibit the same weighting trade-offs as 
experienced in classical control with multiple gain tuning. In order to perform path 
planning, knowledge of the spacecrafts’ working environment is required so that 
coordinate interaction can be determined. With precise knowledge of the working 
environment and each spacecraft’s location, one may be able to specifically command 
and control all spacecrafts' interaction optimally with a centralized control scheme. 
However, one must allow for the exhaustive constraint evaluation, algorithm computation 
time, and universal spacecraft communication requirements which globally optimal 

algorithms require. As the DOF of the spacecrafts’ spatial working environment and the 
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number of spacecrafts increase, the computation time increases dramatically. For 
multiple spacecraft organization, in either 2D or 3D space, there are an infinite number of 
theoretical initial conditions and desired final configurations. This requires calculation of 
the best course of action for each and every spacecraft at a specific point in time and then 
instantaneously commanded them into a precise course of action. In addition, the 
simultaneous communication requirement between the centralized controller and each 
spacecraft is intensive. 

Due to the variable positions and attitudes that a group of satellites may be 
arranged in, a more flexible adaptive control algorithm may be desired. Decentralized, 
adaptive and evolving control algorithms allow for variations in the environment and 
spacecrafts, but are seldom efficient. As the environment changes, the robust controller 
can be designed to function even though the exact scenario was not anticipated. Some of 
these adaptive controllers focus use genetic algorithms to allow each spacecraft to 
function with limited to no global knowledge or communication. These schemes are 
harder to predict and may not even arrive at a solution to a reasonable spatial command. 
As the control system shifts in pursuit of better fitness in multiple spacecraft response, 
there is a delay in commanding. The spacecrafts’ movements often appear clumsy and 
slow. Although with proper development, autonomous and distributed functioning can 
replace strict control and centralization [46]. Simple spacecraft (insect like) self 
organizing theories may be interesting in order to show complex collective behavior. 
However, allowing for emergent behavior of spacecraft is not a risk that most space 
agencies, or companies, are willing to take with their multi-million dollar spacecraft. 
Therefore, it is desirable to have each spacecraft efficiently control their own action 
within the execution of an overarching general command. 

There have been a numerous spacecraft formation control schemes which are 

based on dual level algorithms. The first level determines the overall position and 

orientation of the formation and the second level commands the individual spacecraft to 

achieve and maintain their desired position. The selection of the position and orientation 

of the formation may be determined based on the mission or on the initial position of the 

spacecraft. Both leader-follower and virtual structure formation approaches have 

advantages for rigid formation control [10]. However, neither is effective for multiple 
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spacecraft close proximity operations. Use of a potential function algorithm can 
incorporate the favorable insensitivity to leader perturbation and limited communication 
requirement of leader-follower algorithms with the favorable coordinated behavior of 
virtual structures. The primary attraction point of the potential field acts both as a leader 
and the geometric center of the virtual structure. The attraction point does not have the 
single point of failure issues of a true leader and does not rely on the rigid relative 
positioning of a virtual structure. 

Previous spacecraft work usually limited collision avoidance to predictive states 
which must be include into maneuver planning. These are based on numerically complex 
and computationally intensive uncertainty ellipsoids. The results were collision 
proximity and collision probability indicators that needed to be looped back into mission 
planning decision schemes [47]. This processing may be acceptable for ground stations 
with access to space situational awareness data, but not typically practical for on-board 
control algorithms. The direct inclusion of uncertainty information, such as disturbance 
or sensor covariance, into obstacle’s region of influence eliminates the need for 
uncertainty ellipsoid generation and propagation. Even if these ellipsoids can be 
generated, they are highly dependent on sensor and communication reliability. 

The goal of this research is to develop an efficient control algorithm for use 
during proximity operations for distributed autonomous multiple spacecrafts. The 
distributed nature of the spacecrafts will require that global knowledge is not available to 
each spacecraft [11]. A centralized controller is assumed not to exist, such that each 
spacecraft must perform their portion of the operation with local information and limited 
communications. Each spacecraft must also perform efficient fuel consumption 
maneuvers, while avoiding collisions, and satisfying their commanded goal position and 
orientation. The primary control algorithm is based on APE development from sensory 
inputs. An overview of APE control and related topics follows. 
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B. ARTIFICIAL POTENTIAL FIELD (APF) OVERVIEW 

Artificial potential field functions have been used extensively in robot navigation 
and eontrol [2][3][4][5][6][7][9][48]. The idea of the use of potential funetions for robot 
tasking was pioneered by Khatib [2], Khatib applied the problem to robot navigation for 
obstaele avoidanee with perfeet information. The topologieal ehallenge of the use of the 
potential funetions in feedbaek eontrol laws in order to drive robots was explored by 
Koditsehek [3]. The idea that a eonservative system will settle toward the loeal minima 
of their potential energy, is extended to establish an artifieial potential minima 
eorresponding to an spaeecraft’s goal position. The goal position acts as an attraction 
point for the spaeecraft. Additionally, forbidden positions, sueh as obstaeles, ean be 
assigned suffieiently large potentials values sueh that the spaeeeraft is repelled from these 
regions. The spaeeeraft moves along the negative gradient of the topographie potential 
field to approaeh the minimum artifieial potential value in its environment. Obstaele 
avoidanee with perfeet a priori information was further refined by Rimon and Koditsehek 
[7]. Newman and Hogan [6] extended the applieation of potential funetions for time 
varying goal and obstacles. They showed that APF are effective in simple obstacle 
environments and safer than most path planning algorithms in highly dynamie 
environments. 

For statie environments the artificial potential field and global minimums ean be 
determined off-line. However when the environment is unknown and the presenee of 
obstacles are determined by on-board sensors, the loeal potential field is only known at 
that instant. The real-time motion planning of the spacecrafts can lead to loeal 
minimums. Overeoming the presenee of undesired loeal minimums ean be diffieult due 
to the eomplex shape of potential obstaeles. Although, there are various ways to 
overeome the loeal minimum problem, ineluding defining artifieial potential fields that 
only have a minimum at the goal position and developing methods for eseaping from 
loeal minimum [4]. Due to the limited number of obstaeles in the spaeeeraft environment 
and the dynamics of motion, local minimums do not tend to oeeur. However, as multiple 
spaeeeraft eonverge to a goal position they paek into a spherieal goal region whieh may 
eause eaeh spaeeeraft to ehatter about a loeal minimum. Also ehattering ean oeeur during 
passage between two obstaeles. The repulsion field of the obstacles may influenee the 
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APF based controller to bounce back and forth between the obstacles as it passes between 
them. Precise and smooth motion around obstacles may require an spacecraft to look-up 
pre-computed information about obstacles in a respective region and adjust the potential 
field accordingly [48]. Adjustment of repulsive fields around obstacles and time 
dependent smoothing of potential field was used, along with sliding mode control, by 
Guldner and Utkin [48] in order to navigate robots about fixed obstacles. Refinement of 
the obstacle repulsive field may be incorporated into the control algorithm by using 
minimum communication to periodically update the position of cooperative spacecraft. 

As the spacecraft rendezvous and approach docking, multi-spacecraft control may 
necessitate task assignment among the spacecrafts. Order of docking and servicing may 
matter either due to heterogeneous functions of each spacecraft or efficiency of relative 
motion. The distribution of tasks and computation may require the development of a 
communication protocol. Sliding mode swarm control, with artificial potentials was 
further developed by Gazi [5] [9]. Yang et al. [1] proposed a fuel optimal spacecraft 
formation reconfiguration using multi-spacecraft task assignment. They established a 
communication protocol which allowed for optimal distributed control algorithm to select 
the most fuel efficient task for each spacecraft based on the cost to the entire group. A 
similar communication algorithm seems promising for terminal convergence of 
spacecraft using a potential field based control algorithm. This docking stage with 
multiple spacecraft can also be considered the terminal phase of on-orbit assembly [49]. 
This final stage may require that the local attractive potentials be re-assigned or 
proportionally weighted for each spacecraft. Since assembly sequences can and usually 
are predetermined before launch, it may be possible to design the docking behavior into 
the control algorithm. In addition, predictive and smoothing filters may be applied to the 
potential field as finite control resources allow. 

Artificial potential field guidance was considered for orbital vehicles by C. R. 

Mclnnes in 1993 [50]. It has been expanded to consider distributed control [51], 

autonomous rendezvous with fixed obstacle avoidance [52], autonomous control of on- 

orbit assembly [40], and fuel efficiency constraints for cluster formation [53]. Recent 

application of potential functions in controlling swarms of micro-utility spacecraft also 

shows promise [54][55]. This research explores the use potential function based on 
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velocity error, as opposed to position errors, for controlling a small spacecraft. In 
general, the artificial potential shows promise in stable rendezvous, collision avoidance, 
while maintaining reasonable fuel efficiency. Autonomous and distributed development 
of local potential field algorithm depends on local sensor information. With very little 
communication requirements, cooperative spacecraft obstacles and unknown obstacles 
can be delineated. 

Knowledge about cooperative spacecraft can be used to refine obstacle avoidance 
potential fields. This may be especially useful for terminal stages of the docking. In 
addition, the space environments itself is relatively obstacle free and most potential 
obstacles can be characterized by size. The knowledge that unknown obstacles larger 
than a few meters in diameter are highly unlikely, is useful when estimating an obstacles 
region with only surface detection sensors. In addition, obstacles crossing the orbital 
path will usually be at high enough velocity that collision avoidance maneuvers are not 
necessary or possible for a small satellite with local sensor information. These sorts of 
cosmic collision avoidance require guidance control from higher level system with access 
to additional information. 

The determination of the local potential filed is dependent on the local sensor 
information and limited communication of cooperative spacecraft current states. 
However, the APT control algorithm can be refined to be made more efficient. The 
initial definition of the APT control algorithm is usually forced to choose between 
absolute collision avoidance, stable convergence, or optimal fuel efficiency. Exploration 
and refinement of a time varying control algorithm, based on analytic potential functions, 
should be able to ensure collision-free, Lyapunov stable, and near-optimal fuel efficiency. 
True fuel optimization would require global knowledge, intensive computations, and 
centralized coordination between all group spacecraft. Local knowledge and limited 
communication/coordination with function constraints will not give the fuel minimized 
solution. However, simple iterative optimal schemes may be computationally practical 
while improving efficiency. 

The control of the spacecraft is based on recursively stepping through an APL. 
The potential field is time varying as the sensor readings are updating. This potential 
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field is evaluated by the eontrol algorithm in order to determine an instantaneous step size 
and direetion for the eontrol actuators. The step size is translated into thruster firing 
duration. The orientation of the thrusters must allow for motion along the opposite of the 
potential field’s gradient. Steepest Descent algorithms restrict the search direction 
orthogonal to the contour lines. Refer to Appendix A for a full discussion of optimal 
search algorithms. In some instances convergence performance may be improved if the 
search direction is adjusted to be conjugate. This usually means more steps are required 
for convergence. In fact, convergence is not guaranteed. However, having a known 
location of the goal position will allow us to force convergence by using this position 
information to limit the influence of any possible local minimums. A fuel efficient time- 
varying potential field control algorithm can be developed for autonomous, distributed 
spacecraft close-proximity operations, to include rendezvous and docking maneuvers. 


C. PATH PLANNING CONTROL SCHEME 

Typical path planning algorithms presume a priori knowledge of the spacecraft’s 
geometry and environment. The algorithm establishes an optimal trajectory based on 
some cost function. Then the controller tracks the spacecraft along the trajectory al 
closely as possible. An example of a current satellite trajectory optimization model is 
presented in [45]. These types of optimal controllers are designed to minimize fuel for 
during large spacecraft maneuver or long term formation maintenance. Designing these 
fuel efficient maneuvers for a single vehicle is time consuming, involving iterative check 
to ensure minimization of fuel consumption. This translates into heavy processing 
requirements that seem to be beyond the current on-orbit capability. These optimized 
controllers seldom account for multiple spacecraft and the spatial interactions between 
spacecraft. They typically are only concerned with moderately static formation 
maintenance, such as in [56]. If they do address collision avoidance, by defining a 
relatively larger safe distance constraint in their performance evaluation, the results are 
usually over constrained and one can question weather these methods are the most 
effective in time or energy. For an example of a trajectory planning algorithm related to 
collision avoidance during docking; refer to Breger and How’s research [57]. The 
existence of numerous constraints often results in trajectory planning that is not obvious 
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and computationally intensive. Additionally, in a dynamie environment with 
maneuvering obstaeles or limited environmental knowledge the optimal eontrol must be 
iteratively eomputed and updated. For autonomous distributed eontrol there is a mueh 
simpler solution. 

Instead of foeusing on the eomputational intensive optimal solution for every 
possible multiple spaeeeraft seenario, a feasible, fuel effieient, and safe solution for 
robust, autonomous, and distributed elose proximity operations is desired. In order to 
satisfy these requirements, a suitable potential funetion based closed loop feedbaek 
eontrol approaeh was researehed. This eontrol strategy allows for eomplex multiple 
spaeeeraft behaviors to arise from simple single vehiele aetion. It is generally true that as 
spaeeerafts move from their initial position into the new eommanded position, or 
eonfiguration, they must avoid eollision with eaeh other and potential obstaeles. This is 
espeeially true in orbital operations, where obstaeles may inelude portion of the launeh 
vehiele and unresponsive spaeeeraft. This eollision avoidanee in path planning has 
usually been negleeted for spaeeeraft, due to the large working spaee and relative 
distanee in whieh spaeeeraft usually operate. However, for proximity operations of 
spaeeeraft eonstellations, formations, and swarms/elusters the eontrol of the spaeeeraft 
must also eonsider the relative motion of eaeh spaeeeraft. As seleeted spaeeerafts 
eonverge or disperse in a eontrolled manner, eaeh must maintain a eollision free path. If 
the relative motions are too fast then eontrol forees and torques may reaeh saturation and 
risk overshoot and eollision [49]. If the relative motions are too slow the orbital 
perturbations in elose proximity will eause exeessive eontrol to maintain safe positioning. 
The seleetion of veloeity dampening terms in the potential field is dependent on the 
spaeeeraft’s dynamie environment the aetuators responses [3]. The total effeet of the 
spaeeeraft motion ean be spring-like as the spaeeeraft is pulled toward goals and pushed 
away from obstaeles. 

Proper seleetion of attraction and repulsion potential values for the multiple 

spaeeeraft environment offers promising navigation results. The seleetion of attraetion 

potentials approaehes the fuel and time optimal solutions. The resulting maneuvers tend 

to differ from a true minimum fuel maneuver moderately due to the short time eonstraints 

and eonsideration of multiple parameters. The seleetion of repulsion potentials around 
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obstacles and group spacecraft serves as collision avoidance. Sequential calculations of 
the repulsion potential for converging spacecraft can be considered a near-miss penalty 
function [53]. The superposition of the goal and obstacle potentials function results in a 
local artificial potential field for the spacecraft. By descending the lowest path toward 
the minimum values, the spacecraft approaches the goal similar to a stream flowing down 
a hillside. 

The close proximity operation may be considered complete once the total velocity 
of the spacecrafts falls below a certain threshold within a desired spatial region. The 
potential function is not static, due to its dependence on the relative position of 
spacecrafts and obstacles. Also, sensor uncertainties, dynamic perturbations, and non¬ 
ideal actuators affect the rate of change of the potential function. Therefore, as the 
spacecrafts converge on the final close-proximity goal position they may be subject to 
chattering phenomena [5]. Resolving chattering usually requires rigorous analysis of the 
control commands and model uncertainties. 
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VI. APF CONTROL ALGORITHM APPROACH 


The control of multiple spacecraft during close proximity operations by 
implementing a recursive APF based algorithm was explored. A robust and unique APF 
with collision avoidance control algorithm was developed for close proximity spacecraft 
missions. The APF relates position based potentials to modify the Chase spacecraft 
desired velocity. The collision avoidance used the obstacle positions to dampen the 
Chase spacecraft velocity components in the direction of obstacles. The use of Chase 
spacecraft positioning and velocity relationships in the APF gains allows for 
straightforward adjustment for spacecraft constraints. Additionally, relative positioning 
logic limits the effect of obstacles for precision missions, such as docking. Successful 
simulations for multiple spacecraft during simultaneous maneuvering for convergence, 
rally, rendezvous, docking operations were conducted. Refer to Chapter VIII for 
simulation analysis details. 

The artificial potential field of each spacecraft is determined by the arithmetic 
superposition of the goal and all obstacle potential functions in its working area [6]. The 
APF is an intuitive geometric method of determining the behavior of groups of robots 
[58]. The overall potential field will serve as the performance surface for the control 
algorithm, of the form 

y=y^+K (6.1) 

where is the attractive potential of the goal point and is the repulsive potential of 

obstacles [59]. Selection of the potential functions is critical in ensuring smooth potential 

fields that are stable and provide the desired performance. It may be useful to think of 

APF control algorithms as being geometrically based, with ranges to the goal and from 

obstacles being the primary influences on performance. One strategy to selecting 

potential functions is to base them on the desirable characteristics of Lyapunov functions; 

refer to Chapter VTA. Selection of a potential function with the minimum at the goal 

position can be used to attract spacecraft; refer to Chapter VLB. In addition, potential 

functions with relatively large maximums can be used to push away spacecraft; refer to 

Chapter VI.C. Generally, the most elementary functions that can characterize an obstacle 
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will yield the best results. In 3D spaee, spherieal or ellipsoid shape obstaele fields are 
effective in minimizing local minima. However, the obstacle must have a repulsive field 
in order to ensure collision avoidance. The challenge is to define a region of influence 
around obstacles. If the obstacle region of influence is too large, it will severely limit 
motion and causes excessive avoidance. On the other hand, if the obstacle region of 
influence is too small then collision with the boundary is possible. 

A. LYAPUNOV STABILITY FOR POTENTIAL FUNCTIONS 

The stability of a dynamic system is of interest when engineering a control 
algorithm. The description of potential function stability with respect to a dynamic 
system is generally heuristic. One conceptual method of generating a stable control 
algorithm is to apply the concept of Lyapunov stability. The definitive papers of Kalman 
and Bertram [60] [61] applied the Lyapunov method to a wide variety of control 
problems. The basic concept is that for any isolated system, if the rate of change of the 
energy is negative except for a single minimum equilibrium point, then the energy of the 
system will continue to decrease until the system arrives at the equilibrium point. This 
idea can be expanded into a mathematical definition. 

A dynamical system, with states (x) is asymptotically stable about some 
equilibrium point {x^) if there exists a Lyapunov function ,V(x), with the following 
properties: 

1) Lyapunov function is positive definite, except at equilibrium point. 

y (x) > 0 except at x^ x^ (6.2) 

2) Rate of change of Lyapunov function is negative definite, except at 
equilibrium point: 

y(x)<0 except at x^x^ (6.3) 

3) Lyapunov function is zero at equilibrium point: 

y(xj = 0 (6.4) 

4) Rate of change of Lyapunov function is zero at equilibrium point: 


70 



nxj = 0 (6.5) 

The asymptotic stable equilibrium point can be considered global asymptotically 
stable if it is the only equilibrium point in the domain and the Lyapunov function is 
radially unbounded. A more formal, rigorous, and complete mathematical discussion of 
stability can be found in [62]. 

The Lyapunov function can be extended into a term referred to as a potential 
function. The potential function can be thought of as a topographical representation of 
the state vector of a dynamical system with respect to a desired goal. The system can be 
driven toward equilibrium based on the rate of change of the potential function. A 
sample potential function, shown in Figure 6.1, converges in the direction of the negative 
rate of change. Figure 6.1 shows a 3D plot of quadratic function with the equilibrium 
located at the origin of the R-S plane. This quadratic function meets all of the above 
properties of a Lyapunov function. This type of equilibrium point is often referred to as a 
stable node, or minimum. 



S-axis R-axis 

Figure 6.1 Simple Quadratic Function. 
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The same quadratie funetion is shown on a eontour plot, Figure 6.2, with lines of 
equal potential appearing as eoneentrie ellipses. Vectors of positive slope point toward 
the outside of the contour ellipses and away from the minimum at the origin, while the 
vectors of negative slope point toward the inside of the contour ellipses and in the general 
direction of the minimum point. The steepest slope is characterized by the gradient 
vector, which is perpendicular to each contour ellipse. The negative gradient is the path 
of the Steepest Descent toward the minimum. Refer to Appendix A for discussion of 
Steepest Descent search method. 



Figure 6.2 Simple Contour Plot. 


Since the potential function based on the state vector is a Lyapunov function, it 
maintains the same characteristics listed in equations (6.2) - (6.5). For instance, any path 
following the negative rate of change will arrive at the equilibrium point. 
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B. GOAL POTENTIAL FUNCTION 


The goal potential is the minimum potential in the spaceeraft working area. This 
equilibrium point serves as the global attraetion point for the Chase spaeeeraft. A simple 
quadratie goal potential funetion based on the relative position of the Chase spaeeeraft to 
the goal position is 

( 6 . 6 ) 


where is the relative position veetor of the Chase spaeeeraft, is the relative position 
veetor of the goal, and is a non-negative sealing parameter. The quadratie nature of 

this potential funetion V ensures that it is positive semi-definite, such that it is positive 
for all values of ^ and zero only at ^ = fg • This potential function takes advantage 

of the useful characteristics of both performance functions and Lyapunov functions. 

Next, ensure that this goal potential function rate of change is negative semi- 
definite. Differentiating the goal potential function, equation (6.6), yields the following 


where v is the velocity vector of the Chase spacecraft. The control algorithm ensures 
that the rate of change of the potential function is negative for any ^ ^ and zero when 

the Chase spacecraft arrives at the goal position, r^=f^. This requires that Chase 

spacecraft velocity is maintained along vectors toward the negative slope of the potential 
function. The desired velocity of the Chase spacecraft may be determined as 



( 6 . 8 ) 


where is a positive parameter or a magnitude shaping function [40] and W is the 
gradient of the potential function. The normalized gradient function points the Chase 
spacecraft in the Steepest Descent direction, so that V is always negative semi-definite. 
The selection of determines the convergence of the control algorithm by relating the 
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potential function to a convergence velocity. Large values of cause the algorithm to 

converge quickly toward the area of the goal position, but oscillate around the actual goal 
position. Small values of ensure slow steady convergence toward the goal position in 

a damped manner. This is the more desirable of the possible behaviors for spacecraft 
rendezvous; however convergence to the goal may take more iterations and longer 
amounts of time. The maneuver is controlled along the Chase spacecraft’s relative 
velocity vector coinciding with the Steepest Descent approach. 

The shaping function, k^ , can be directly related to the actuator, such as on-off jet 

thruster performance. The shaping function can be tuned for different spacecraft by 
taking into consideration the minimum thrust pulse duration. The thruster mapping 
algorithm for transforming the control algorithm commands into the thrust components is 
not considered as part of the APT development. This mapping is dependent on the 
thruster configuration and performance on the specific spacecraft. 

C. OBSTACLE POTENTIAL FUNCTION 

An obstacle potential aids the spacecraft to avoid collisions. The obstacle 
potential is a useful tool for avoiding special regions. For instance, free orbiting objects 
and regularly maneuvering spacecraft may be avoided. However, collision avoidance is 
not intended for active evasion of high-speed celestial objects. The location of an 
obstacle may be from a priori information, such as known space debris in a particular 
orbit. Otherwise, obstacle location information may come from onboard proximity 
sensors. If the obstacle is known to be another spacecraft then its size and center of mass 
can be estimated in order to refine the obstacles potential function. The center of mass 
and geometric center of an obstacle can usually be assumed to be the same. Otherwise, a 
Chase spacecraft’s local range sensor may be the only information on which to base the 
obstacle potential function. In this case, the obstacle potential function region of 
influence must be larger or equal to the actual physical region occupied by the obstacle. 
It is assumed that the slight loss of fuel efficiency due to over estimating the size of an 
obstacle is worth the effective execution of collision free navigation. 
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The goal potential and attraetive potential ean be formulated separately, but must 
have properly sealed parameters in order to allow for the desired system motion [58]. An 
obstruetion potential function, based on [52] and [40], is 

(r -r (r -r ) 

^ c o' ^ c o' 

y, = 4 e " (6.9) 

where is position vector of the obstacle, is a scaling parameter for the size of the 
potential function, and <j is the standard deviation of the region of influence. This 
potential function is a Gaussian with and <j selected to ensure that the obstacles 
region of influence (T>„) is equal to or larger then the actual dimensions of the object to 

be avoided. Numerous candidates could be selected for the obstacle avoidance function, 
such as spherical power-law and super quadratic functions. However, these functions are 
more complex to define and generate then the Gaussian function and do not guarantee 
better performance [40]. This complexity would require more a priori of the obstacles, 
which is not assumed in this research. 

Previous obstruction potential functions have been scaled in order to be equal or 
greater than the goal potential at the initial position [52]. This selection is 
mathematically useful so that the obstacle potential is large enough to ensure that the 
Chase vehicle avoids the obstacle region. However, a potential function based only on 
relative positions may result in circular local minimums around obstacles. Also in the 
presence of multiple obstacles, the global minimum may be shifted due to the 
superposition of obstacle potentials [52]. If the APT function is purely position based 
then velocity variations may cause unacceptable overshooting position oscillations in the 
region of obstacles and the goal location. Relating the potential function with Chase 
vehicle’s desired velocity may eliminate such performance limitations. The proper 
consideration of velocity in the application of APT methods is essential for application on 
practical systems with limited actuation. The velocity relationship can be included of as 
a velocity error potential function, as in [54] and [55], or as coupled relative state APT, as 
in [63] [64]. The coupled relative state APT is developed in this research as a position 
based APT which drives the desired velocity. The magnitude and slopes of the potentials 
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are driven by the relative position between the Chase spacecraft and the goal and 
obstacles. A novel selection of the APF relationships and parameters for the multiple 
spacecraft was determined. 

D. APF FOR MULTIPLE SPACECRAFT PROXIMITY OPERATIONS 

The APF can be refined for application to the multiple spacecraft close proximity 
operation. The goal and obstacle APF functions for relative spacecraft dynamics was 
selected for precision of relative position, steady convergence, and control efficiency. 
The simplified APF control block diagram is shown in Figure 6.3. 


Desired (r, v) 



Figure 6.3 APF Control Block Diagram. 


The goal, or attractive, potential remains the same as quadratic equation (6.6), 
with the positive scaling parameter = 1 / ~ ^ ^ that 


V =-r 

g 2 


( 6 . 10 ) 


This goal potential function results in a linear scaling of the range. Range maintains a 
positive semi-definite relationship based on the vector normalization. 
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Application of APF to the spacecraft environment requires very careful 
considerations of acceptable relative velocity. Stable relative orbital dynamics and 
limited spacecraft translational actuation seriously limits acceptable velocity magnitudes 
and variations. Due to these spacecraft relative dynamics and control relationship, the 
goal APF is used to determine the spacecraft desired velocity in an exponential manner. 
This acts as a damping of the rendezvous dynamics necessary for position precision in the 
presence of actuator limitations. This goal potential is used to scale the desired relative 
velocity of the Chase spacecraft, as in equation (6.8). The desired velocity magnitude 
shaping function is a positive semi-definite shaping function. The shaping function, 

kg , was determined to be 
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^ r. . ^ 

mit 
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V 'max J 





( 6 . 11 ) 


where the ratio of the initial range, , and maximum range, , linearly scales the 

maximum relative velocity, . The variable, bg , is used to shape the exponential 

decay. The value allows for slower decay, and can be scalable to the initial conditions of 
the Chase spacecraft, such as 


= 




v"',/ 


V ^init ) 


( 6 . 12 ) 


where is a positive scaling constant. In this research dg=\l was selected for 

favorable velocity decay characteristics. The initial position determines the initial desire 
velocity magnitude which is exponential shaped as the Chase spacecraft approaches the 
goal position. The resulting Chase spacecraft’s desired velocity, based on the attraction 
toward the goal position, is 
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(6.13) 
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This desired veloeity shaping funetion deereases the Chase spaeeeraft desired 
veloeity as it approaehes the goal. Based on an assumption of zero starting relative 
veloeity, the initial veloeity transient is often large and eauses the eontrol aetuator to 
saturate. In order to avoid this saturation an exponential ramping funetion, , ean be 
added to Equation (6.13), sueh as 



\ 'max J 


V„ 




(6.14) 


where this veloeity ramping term ean be represented as 


k 


R 



(6.15) 


For this researeh the veloeity ramping eonstant, was seleeted to be =10. This 

allowed for more eomparable performanee with the eontrol algorithm developed in 
Chapter VII. This ramping term only influenees the initial veloeity transient by allowing 
a more gradual inerease related to the start-up of the eontrol algorithm and the initial 
range from the goal. For instanee, maneuver of approximately 100 meters result in a 
ramp up to the maximum veloeity is about 60 seeonds. 

The aetual relative veloeity is subtraeted from the desire veloeity to determine the 
Av required by the eontrol effort. This desired ehange in veloeity is used to determine 
the Chase spaeeeraft’s desired eontrol aetuation, in terms of aeeeleration 



The goal potential allows for eonvergenee to the goal position, however an 
obstaele potential is required to avoid eollision with other spaeeeraft and sensed objeets. 
The obstaele avoidanee funetion is best thought of as a damping of the geometrie forees 
toward an obstaele, sueh as relative veloeity and aeeeleration. Due to the manifold of 
relative dynamies, aetual pushing away from an objeet ean lead to instability. The 
obstaele potential seleeted to be a Gaussian funetion of the form 
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(6.17) 


where , is the distanee from the Chase spaeeeraft to the eenter of the obstaele, is 
the region of influenee from the eenter of the obstaele, cr is the standard deviation of the 
bell-shaped eurve, and is a positive funetion whieh serves as a repulsive potential 

magnitude shaping function. The repulsion potential curve is a smooth function that 
increases from the boundary of the region of influence to the surface of the obstacle. A 
smooth transition is ensured by subtracting the value of the Gaussian function at the edge 
of the region of influence. The employment of a Gaussian function allows for direct 
implementation of uncertainty in an obstacle state. The shaping parameter is used to 

ensure that the value of the repulsive potential at the obstacle surface is equal to the initial 
attraction potential. The obstacle shaping parameter used in this research is 
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(6.18) 


where is the exterior surface radius of the obstacle, such as another spacecraft. The 
region of influence of the obstacle, , and the standard deviation, cr, are functions of 
the size of the obstacle, . This selection of A^ ensures that the value of equals the 
initial value of at the surface of the obstacle. The region of influence, d^, is 

determined from the size of the obstacle, the velocity of the Chase spacecraft, v , and the 
maximum acceleration, , allowed by the control actuation. The minimum stopping 

distance, 0^,^^ required by a spacecraft is determined as 


D 


Stop 


V 




(6.19) 


This functional link allows the speed and responsiveness of the Chase spacecraft in the 
rendezvous region to determine the buffer distance which is required for obstacles. 
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( 6 . 20 ) 


D ^ d (l + D , ) 

o o y o stop j 

where the eonstant multiple faetor, allows for a smooth breaking region for avoiding 
impaet with obstaeles. For this researeh, d^=3 was seleeted for favorable avoidanee 
eharaeteristies. The standard deviation cr is seleeted so that the obstaele surfaee is 
within one standard deviation as the spaeeeraft veloeity approaehes zero, sueh that 



This relationship allows a reasonable safety region around obstaeles and a smooth 
Gaussian repulsive potential funetion. The veloeity ehange due to applying maximum 
thruster aeeeleration is shown in Figure 6.4. From an initial veloeity and a maximum 
aeeeleration the minimum stopping distanee is determined by equation (6.19) and the 
resulting veloeity during the braking is shown as a solid line. The aetual braking ean not 
be steeper than this maximum eondition. Using the obstaele shaping parameter, equation 
(6.18), and the standard deviation, equation (6.21), a smooth Gaussian veloeity ehange 
ean be eommanded, as shown as a dashed line. Notiee that both funetions are equal to 
zero at the obstaele’s outer surfaee. For praetieal applieation the Gaussian must not be 
steeper than the maximum braking eondition at any point. Also, the area of the 
maximum stopping funetion must be less to the area under the Gaussian eurve. This is 
aehieved in the algorithm by using the three standard deviations with an additive term. 
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Figure 6.4 Spacecraft Velocity in Obstacle Region. 


However, application of this position based obstacle repulsion potential requires that it is 
related to the attractive potential. This obstacle potential is used to modify the desired 
relative velocity of the Chase spacecraft, as in equation (6.8). The desired velocity 
shaping function is a positive semi-defmite shaping function used push the Chase 

spacecraft away from the obstacle. The shaping function, , was determined to be 


k V 

k=^ 


imt 


( 6 . 22 ) 


This unique shaping function allows for k^ being equal to k^ at the surface of the 

obstacle. The resulting Chase spacecraft’s desired velocity due to the obstacle, based on 
the repulsion away from the obstacles position, is 
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(6.23) 
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The obstacle shaping function dampens the Chase spacecraft approach in the 
direction of the obstacle as spacecraft advances toward the goal position. Proper 
selection of damping based on components of the relative velocity vector results in a 
smooth flow around spherical obstacles. Additionally, this insight and use of damping 
shaping functions should eliminate of any possible narrow passageway oscillations, such 
as discussed by Masoud [65]. The attraction velocity vector of the goal is toward the 
goal position and the repulsive velocity vector of obstacles away from each obstacle. The 
effect of each obstacle within the spacecraft’s spatial region is summed in order to 
achieve the total repulsive force on the Chase spacecraft. The total control force is 
determined by vector addition of the potential derived velocity minus the current actual 
velocity vector of the Chase spacecraft, as 
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^total 
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(6.24) 


where the number of obstacles is not limited, but can result in Chase spacecraft not being 
able to arrive at the goal location. Obstacles may be other spacecraft or stationary 
exclusion zones. The other spacecraft are generally moving relative to the Target and 
Chase spacecraft. In this research, other spacecraft are typically additional Chase 
spacecraft converging toward a goal within the same region. Stationary obstacles are in a 
fixed position relative to the goal location. These may be extensions of Target spacecraft 
or exclusion zones due to physical objects, such as solar panels, or nonphysical items, 
such as thrust plume or radiation beam regions [52]. 

Selection of the magnitude of the repulsion shaping function must be related to 

the attraction shaping function in order to achieve desired critically damped performance. 

Proper selection of a repulsion shaping function based on the attraction shaping function 

allows for safety in selecting goal positions and efficiency when avoiding obstacles. For 

instance if the region of influence of the obstacle is too small and the slope of the 
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repulsive potential shaping funetion is too steep then a thrust limited aetuator may not be 
able to avoid eollision with the obstaele. On the other hand, if the obstaele region is too 
large then the Chase spaeeeraft may be less effieient in both eontrol effort and maneuver 
duration as it avoids obstaeles. 

An obstaele’s repulsive region of influenee eauses a potential minimum or saddle 
point to oeeur in the area between the obstaele outer region of influenee and the surfaee 
of the obstaele. The loeation of this loeal minimum depends on the obstaeles loeation 
with respeet to the goal position. This loeal minimum ean eause diffieulty if the overall 
potential funetion is the only driving funetion for determining control effort. However, 
the attractive and repulsive desired velocity shaping functions, and respectively, 

allow for velocity damping around regions of concern. This ensured that the Chase 
spacecraft slows as it approaches the goal position and avoids obstacles. These shaping 
functions are also selected so that the desired velocity determined from the repulsive 
potential balances with the attractive potential at the surface of the obstacle. This allows 
for the goal position to be placed in the center of a spacecraft and the control algorithm to 
converge to the surface of the Target spacecraft. This is a vital capability for docking 
algorithms. The control algorithm is precise within a millimeter of the goal position or 
the Targets spacecraft’s outer surface, assuming ideal sensor measurements. 

As multiple spacecraft and obstacles occupy the Chase spacecraft’s region, some 
simple control logic must be applied. First, Chase spacecraft are only influence by 
obstacles when they are within the obstacles region of influence, . Second, only 

obstacles which are equal distance or closer to the goal position than the Chase spacecraft 
are allowed to influence the Chase spacecraft. For instance, the spacecraft is looking 
toward the goal like an automobile on the road. The next spacecraft converging into the 
goal region will then avoid contact with Chase spacecraft. In most cases, other spacecraft 
are simply treated the same as obstacles. However, additional logic is needed if multiple 
spacecraft are converging within the same goal position. The third control logic 
condition uses a safety function, k ^, to modify the desired repulsive velocity between 

maneuvering spacecraft as they approach the goal. This safety influence between 
multiple converging spacecraft ensures collision avoidance while achieving the closest 
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possible convergence to the goal. This safety function is a trade-off between collision 
avoidance and docking. The safety function, , between converging spacecraft 
multiplies the current and results in the following modification to (6.23) as 

=KK— (6.25) 

^co 

where is usually greater than or equal to one. If k^=\, there is an ideal balance 
between and at the surface of the goal location boundary. This may allow contact 

between Chase spacecraft converging at the goal location, which obviously can be 
undesirable. A value of =1.01 is large enough to ensure that multiple spacecraft 
converging upon the exact same goal position do not collide. However, the multiple 
spacecraft rendezvous to the exact same goal position results in a staggered convergence. 
The first Chase spacecraft to arrive converges to the goal position. The next Chase 
spacecraft converges to within millimeters of the first Chase spacecraft. The third Chase 
spacecraft has the additive repulsion of the first two spacecraft and converges to a radial 
position further away. Any additional spacecraft will converge to a safe point slightly 
further away. This staggered cluster may be a desirable result for spacecraft 
rendezvousing to an unknown formation, where additional command maneuvering may 
need to occur. 

For multiple spacecraft docking maneuvers, the staggered cluster effect of the 
additive repulsion may not be desired. In this case, the goal location is an actual Target 
spacecraft. To allow the later arriving spacecraft to converge toward docking, while 
avoiding collision, the safety function, k ^, is adapted to be a decaying exponential of the 
attractive potential based on the goal position, such as 

-(v 

K=l-e ^ (6.26) 

This results in the repulsion due to other spacecraft decaying toward zero as the Chase 
spacecraft reaches the outer bound of the Target spacecraft. In this case, the multiple 
spacecraft converge relatively tightly around the Target spacecraft. Limitations in the 
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Target spacecraft’s outer boundary surface area and local minimums due to saddle points 
may cause some delays for spacecraft which arrive late. This is usually only an issue for 
a second wave of spacecraft which arrive as the first few converging spacecraft are 
settling into position. It is envisioned that each spacecraft is commanded to a specific 
docking port, therefore this clustered convergence is not an operational issue. Also, more 
spacecraft will not typically be commanded to converge toward docking on a spacecraft 
which can not support additional spacecraft. 

Some simple logic is used to ensure that obstacle influence and goal convergence 
is reasonable. For instance, any obstacle outside the relative range of the goal position is 
not allowed to influence the Chase spacecraft. This ensures that late arriving spacecraft 
do not cause docking spacecraft to bump into the Target spacecraft. The Chase 
spacecraft’s collision avoidance motion is dependent on three primary logical conditions, 
as follows: 

1. Chase spacecraft must be within region of influence of an obstacle, such that 
r < D . 

CO o 

2. Chase spacecraft range to its goal must be greater than the obstacle’s range to its 

goal, such that . This allows the Chase spacecraft to be influenced 

only by obstacles closer to their goal, so that a closer Chase spacecraft is not 
disrupted from its goal due to a farther converging spacecraft. The inclusion of 
the term acts as a safety margin ensuring that the Chase spacecraft does not 
clip the side of obstacles as it resumes it free-space motion. 

3. Chase spacecraft range to its goal must be greater than the distance to the 

obstacle, such that -L^/2). Obstacles on the far side of a goal location 

are not allowed to influence the convergence toward the goal. The inclusion of 
the term acts as a safety margin ensuring that other docked spacecraft are 
avoided while having a limited influence. 

These logical conditions limit the collision avoidance in obstacle dense 
environments. The most sensitive of these environments occurs when multiple spacecraft 
are simultaneously converging toward a common Target spacecraft. If the desire is for 
multiple spacecraft to simultaneously dock, then the obstacle influence of other Chase 
spacecraft is decayed near the goal location, as shown in Equation (6.26). The third 
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logical condition allows for the Target spaeeeraft repulsive influenee to deeay along the 
doeking port’s axis. This deeay results in eardioids shapes similar to that diseussed by 
Lopez and Meinnes [52], without the need for additional transformations. The Target 
spaeeeraft’s region of influenee is due to the spherieal, Gaussian based, region of 
influenee with a roughly eonie shaped wedge eut from around the doeking port. 
Although, a full 3D representation is not neeessary sinee an obstaele’s region of influenee 
is primarily range dependent. 

This multiple spaeeeraft APT eontrol algorithm appears to be directly applicable 
to any spaeeeraft in a full range of elose proximity maneuvers. The parameters for the 
potentials and shaping funetions are not system speeifie, so intensive tuning is not 
required. The performanee appears to be robust in the full range of elose proximity 
operations. Maneuver ranges of over one kilometer, doeking preeision within 
millimeters, and relatively high allowable speeds establish improvements in timeliness 
and robustness. Direet eomparisons with previous APT eontrol algorithms for spaeeeraft 
is limited, sinee few whieh have been fully developed and simulated. For initial 
eonsideration refer to previous researeh eondueted by MeQuade [40] or the system 
speeifie algorithm reeently developed by Neubauer [54][55]. 


86 



VII. LQR/APF CONTROL ALGORITHM APPROACH 


The LQR approach can be recursively applied to the multiple spacecraft close 
proximity control. The inclusion of dynamics in the LQR computations allows for 
optimality considerations. Also, there is a desire to combine the APT repulsive collision 
avoidance capability with the refined LQR for spacecraft rendezvous. The developed 
LQR with APF-based collision avoidance allows for efficiency based on the 
environmental dynamics combined with gain weighting and logic to allow for primary 
mission achievement. The simplified LQR/APF control block diagram is shown in 
Figure 7.1. 


Desired (r, v) 



Figure 7.1 


LQR/APF Control Block Diagram. 


A. GENERAL LQR 

The LQR uses the state space dynamics of the system to determine the optimal 
control effort based on solving Riccati equation for a selected cost function with full state 
feedback. A quadratic cost function can take into account control effort for a linear 
system is referred to a LQR. The quadratic cost function is of the general form 
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Jlqr = —I Qx + if Ru+if Nx^dt 


The gain matrixes for the states, Q , and eontrol effort, R , are of primary interest. The 
eoupling gain matrix, N , is typieally zero. These gains are used to solve the Riceati 


equation. 




where S is the solution of the algebraie Rieeati equation. The state feedbaek to 
minimize the quadratie eost funetion ean be determined from this solution. The optimal 
state feedbaek gain, K^qj^ , for a partieular eost funetion is determined by 


S + N^) 


This optimal feedbaek gain is then used to determine the optimal eontrol effort, such as 

u =-R^\b'^ S + N'^)x=-K^q^x (7.4) 

The computation for a limited state LQR can be relatively efficient, as compared to more 
complicated optimal routines. Increasing the magnitude of the state weighting matrix, 
Q , will result in a faster convergence to goal. Increasing the magnitude of the control 
effort weighting matrix, R , will result in more efficiency. The trade-off is between the 
time of convergence and the control efficiency. The gain is mostly a constant gain, 

which decreases as the states converge to the final value. 

The determination of the gain matrices is usually based on the normalizing the 
quadratics along the diagonal. First attempts typically use the maximum allowable 
values of the states and control effort as the initial values, such as 


(^Imax)' 


(■^«max ) 
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with adjusted based on the state performance. Similarly, thei? weighting matrix is 
selected based on the maximum allowable control effort, such as 

A, ^ 

(^Imax ) 

0 

with adjusted based on the control effort demanded. The selection of these gain 

variables require insight into the desired response of the system with modifications based 
on simulation results. 

B. LQR/APF FOR MULTIPLE SPACECRAFT PROXIMITY OPERATIONS 

The LQR control response for the multiple spacecraft maneuver requires 
consideration of system response and control efficiency. System response includes the 
concepts of maneuver duration and precision upon reaching the goal. For the multiple 
spacecraft rendezvous problem, a critically damped relative position response with 
limited control effort is desired. As with all spacecraft maneuvers, control efficiency 
during multiple spacecraft close proximity operations must be considered. However, the 
convergence maneuver is assumed to be operationally significant and must be performed 
in finite time duration. For this research, close proximity maneuver durations of one 
quarter of an orbital period, approximately 30 minutes, are desired. This duration is 
based on the spacecraft starting from an initial relative position of within one kilometer. 
The close proximity maneuver is considered successful once the spacecraft converges 
within a spherical region from its desired goal position. This precision may be much 
greater than the typical one meter used to evaluate most rendezvous maneuvers. The 
intent is to be able to use the developed control algorithm for docking maneuvers. 

The balancing factor between spacecraft relative position and control effort 
efficiency is the relative convergence rate. However, the relative spacecraft dynamics 
causes rendezvous challenges if the relative convergence rate is too slow or rapid. If the 
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rate of convergence is slow the goal position is spirally orbited as the minimal control 
actuation is used. The slow converges can dramatically increase the maneuver duration 
as the spacecraft approaches close to the goal position. The long duration within close 
proximity of multiple spacecraft can unnecessarily increase the danger of collision due to 
perturbations. On the other hand, if the rate of convergence is too rapid the limited 
actuation will result in a collision danger due to relative position overshoot. Even if 
collision is avoided initially, the spacecraft continues to orbit through the goal position 
due to overshooting oscillation effects, often referred to as a rubberband effect. This 
limits both the close proximity maneuver precision and efficiency, with the resulting 
motion of the Chase spacecraft oscillating about the goal position. 

Taking the above issues into consideration, the attractive potential of the APE 
algorithm is replaced with an iterative EQR for multiple spacecraft in proximity 
operations. Collision avoidance was incorporated based on the APE-type shaping 
function. The EQR becomes the driving control toward the goal and the APE based 
collision avoidance is successfully applied. The resulting EQR with APE-based collision 
avoidance incorporates linearized dynamics for free-space optimal convergence and 
utilizes avoidance based on obstacle geometric relationships. The developed multiple 
spacecraft EQR/APE control algorithm is control effort efficient and effectively avoids 
collisions while successfully conducting a wide range of close proximity operations. 

1. LQR Attractive Component 

The iterative EQR gain matrices for multiple spacecraft were selected after 
evaluating the APE control responses. The gains were determined for six states, position 
and velocity along each axis, along with control effort along each axis. The resulting 
gains matrixes, based on equations (7.5) and (7.6), are as follows. 
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The relative position error along each axis is equally weighted, so that 


X — \ —7 — r 

'^max -^max "^max eg 


(7.9) 


The selection of maximum distance as the current distance to the goal, r , allows position 
to become more important as the spacecraft approaches the goal. The relative velocity 
error along each axis is also equally weighted, such that 


X = V —7 — 

'^max y max "^max 


^r. ^ 

imt 


V ^ } 

V max / 


(7.10) 


91 



This desired veloeity is determined by sealing the maximum relative veloeity, by the 
ratio of the initial range, , and maximum range, . For this researeh, the maximum 
relative veloeity was conservatively selected as 

This selection of velocity limits the transients due to the initially neutral relative velocity 
and also limits the convergence rate for safe operations. However, this is 

approximately ten times greater than the relative velocities allowed by previous APF 
proximity control [40] [54], 

With the denominators of the diagonal gains selected for desired maximum 
position and velocity. The numerator position error and velocity errors along each axis 
are equally weighted. Therefore, the gains are all equal and scaled as the Chase’s 

distance to the goal position converges 

^Gi ~^Q2 ~^Q3 ~^Q4 ~^Q(, “Cg (’7-12) 


The actuator control effort is the acceleration imparted due to the translational 
thrusters. For this research, the thrust along each axis is limited to a maximum 
acceleration of 


=^^ = 0.01 mls^ 
m„ 


(7.13) 


based on a thrust force of 1.0 Newton and a spacecraft mass of 100 kg. This maximum 
acceleration is the maximum control effort along each axis. 


u =u =u = a 

xmax ymax zmax max 


(7.14) 


The scaling of the control effort is also scaled as the spacecraft relative position 
changes. Therefore, the numerators of the control effort matrix diagonal gains are 


A, =/^R. =PR, =Cg 


(7.15) 
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A minimum scaling factor for the numerator ean be selected such, as the range to 
goal approaehes zero, that numerieal problems and ehattering is avoided. For instance as 
the approaehes zero the value of may be limited to some nominal minimum value, 

sueh as r > 0.05 m . 

^6 


The control effort resulting from the LQR algorithm, refer to Equation (7.4), ean 
be thought of in terms of a desired aceeleration, sueh as 


‘^LQR 


^LQR ^ 


(7.16) 


where v is a generalized representation of the Chase spaeeeraft’s position and veloeity 
states. 


2. APF-Based Collision Avoidance Component 

The LQR drives the eontrol effort based on the relative spacecraft systems 
linearized state dynamies. The LQR algorithm eontrol effort varies the Chase 
spaceeraft’s position and veloeity states in a more complieated manner then the previous 
geometrie based APL. This more complicated relationship requires a modifieation to 
both veloeity and acceleration in the region of influence of obstaeles. The component of 
the Chase’s veloeity in the direetion of the obstaele is determined as 


= 


^ogset ) f 


r 

CO 


K^coJ 


(7.17) 


where is the veloeity offset eorreetion due to the faet that relative stationary 

obstacles aetually maintain a veloeity in the orbital plane. This eorreetion is eritieal for 
potential sub centimeter preeision maneuver, sueh as doeking. Similarly, the eomponent 
of the acceleration in the direetion of the obstacle is determined as 


= 


rco^^LQR 


f- \ 


K^coj 


(7.18) 
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Using the vector components of the velocity and acceleration in the direction of obstacles 
helps limit superposition issues. The resulting iterative spacecraft control algorithm is 
driven by optimal LQR cost convergence, with associated dynamics, and smooth 
collision avoidance responses. 

The APT obstacle potential parameters, represented in (6.22), can be combined to 
generate a Gaussian function which is unity at the obstacle boundary. This function 
becomes the LQR velocity shaping function due to obstacle position. 


k 

V 



(7.19) 


This gain, , is multiplied by the negative relative velocity to ensure the Chase 
spacecraft slows to zero at the boundary of the obstacle. 

Next, the attractive acceleration due to the LQR recursive function is shaped. 
There is no change to it when the Chase spacecraft is outside obstacle regions of 
influence. However, if the Chase is within the region of influence then acceleration 
toward the obstacle must be decreased. The LQR acceleration shaping parameter is 
selected as 


k 


a 



(7.20) 


where the positive constant, , is used to establish the parameters rate of decay. In this 

research = 1 was selected in order for the decay directly related to relative position. 
The parameter is multiplied by the negative component of LQR acceleration to ensure 

that the LQR derived control effort does not drive into an obstacle. Finally, the safety 
shaping parameter, from Equation (6.26), is modified to replace the potential function 
with the Chase spacecraft’s range from the goal. 


k 


S 



(7.21) 


94 



The safety function allows the obstacle repulsion to decay faster as the Chase 
spacecraft approaches the goal position. This function enables precision maneuvers, such 
as docking, in regions where the relative ranges between spacecraft are small. If the 
obstacle is the Target spacecraft then the safety function ensures that the Chase only 
approaches in the vicinity of the docking port. 

The overall control effort for the multiple spacecraft LQR with collision 
avoidance is 


n 

« = ^LQR - S ((^v K ^co ) (7.22) 

obs=0 

The control algorithm only decreases velocity and acceleration toward obstacles. 
It does not actually push away from obstacles. This is useful in maintaining the relative 
stability in a bounded system. Fortunately the relative dynamics result in forces which 
help the control algorithm escape local minima in densely packed obstacle regions. The 
consequence is similar to that achieved by APF wall-following methods [4]. The 
efficiency gained by the LQR derived control effort is more significant when 
implemented in a limited number of obstacles environment. 


C. SENSOR NOISE AND MODEL UNCERTAINTY 

In actual applications, the potential field and related gradient are estimates 
dependent on noisy sensor measurements and spacecraft modeling uncertainties. A 
realistic system model typically has imperfect initial estimates and dynamic disturbances. 
Measurements of the system states are limited and imperfect. These issues can be 
considered in the simulation by applying variance as a measure of the uncertainty. The 
current state of the spacecraft and the latest sensor measurements can be filtered in order 
to determine the best estimate of the spacecraft. For instance, the Kalman filter 
determines this estimate in a minimum mean square error manner [66]. The Kalman 
filter can be implemented as predictor-corrector algorithm. The filter predicts the next 
position of the spacecraft and the value of the next sensor measurement. Next, the next 
sensor measurement is used to determine the residual error between the predicted and 
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actual position. This residual is multiplied by a weighting gain and used to eorreet the 
predieted state to determine the eurrent estimate. 

The number, type, and loeation of sensors are dependent on the degree of 
maneuver preeision required and knowledge of aetuation responses. For simple doeking 
meehanisms sueh as that used on Orbital Express operations arehiteeture, the eonneetion 
of spaeeeraft to a hardened doeking attaehment ring does not require eomplex 
manipulator eoordination. The entire spaeeeraft is positioned and oriented for doeking. 
The need for additional sensors is limited to the end-point of the simple doeking 
meehanism. The aetuator eontrol is usually limited to translational thruster and attitude 
reaetion wheels (or possibly jets). Thruster aetuation variations, sueh as those mentioned 
in Chapter III.E.5., ean be due to thermal ehanges and fuel supply levels. 

The use of manipulator arms may be neeessary for some doeking and assembly 
seenarios. This is espeeially important for eomplex on-orbit assembly of struetures, sueh 
as solar arrays and power stations. As the assembled structure inereases in size the 
flexibility, thermal, and other disturbanees aeting on the body assembly requires more 
detailed analysis. If the doeking and assembly is eondueted with free flying manipulator 
space robots, then the modal struetural oharaeteristies of the struetures must be aeeounted 
for with high preeision and minimum sensors. Researeh at MIT [67] shows promising 
progress for extension of the Base Sensor Control (BSC) method for estimation of 
aetuation forees and torques on spaee robot manipulators with limited sensing. Both the 
foree and torque on a seven DOE bi-arm manipulator ean be effeetively measured with 
only one six-axis foree/torque sensor [67]. The savings in sensor hardware eomplexity is 
highly desirable. 

Even doeking/assembly missions with manipulator arms must get within range for 
manipulator arms to operate. These elose proximity operations must be eondueted within 
the safe operating ranges established for small spaeeeraft. The sensors must be preeise 
enough to measure state eharaeteristies at least to the level of preeision required for the 
spaeeeraft operation. In addition, the aetuators must also be able to operate in a manner 
whieh allows for the desired operational preeision. This may also restriet whieh aetuators 
ean be used at eertain distanees and orientations. Eor instanee, in elose proximity to 
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another spacecraft thruster firings should be restricted to avoid plum impingement on the 
sensors and solar panels [68]. These restrictions are dynamic, since they are dependent 
on the relative position and orientation of the spacecraft. As with any parametric used as 
a constraint or optimization cost function, trade-offs in performance must be properly 
weighted. For instance, selection of different collision avoidance, fuel efficiency, plume 
impingement and time duration considerations will vary the trajectory followed by the 
spacecraft. Each of these metrics may have different levels of sensitivity to model 
uncertainties and sensor noise. These considerations in trajectory and path planning raise 
a large level of concern with the any multiple parameter cost function’s claims on 
optimality. To avoid this level of ambiguity, the close proximity LQR/APF control 
algorithm uses only relative spatial considerations to determine the direction and rates of 
spacecraft motion. 

Estimations of motion between any object, target or obstacles can be improved by 
associating sensor measurements with the actual object position. Adopting the general 
convention that the target states information is assumed to be known, these relative 
motions can be determined by simple vector arithmetic. The difference between this 
expected and measured sensor data is considered the innovation [69]. The covariance of 
the innovation can be directly utilized in determining the collision avoidance thresholds 
for obstacles. These uncertainties are immediately applicable to the Gaussian based 
collision avoidance functions and the EQG algorithm which can incorporate process and 
measurement noise as Gaussian white noises with covariance. The use of Gaussian based 
collision avoidance lends itself to direct inclusion of any measurement noise in to the 
obstacle relative range and velocity. Therefore, the overall control algorithm structure 
allows for convenient inclusion of known, or estimated, sensor uncertainties. 

D. STABILITY AND ROBUSTNESS 

Stability and robustness of collision avoidance control algorithms are not 

guaranteed. The stability is highly dependent on the reference frame and specific 

conditions applied. Generally speaking, once in orbit all passive objects are stable within 

the gravitational region of the primary body. Only by adding energy can the satellite 

escape the primary body. However, once the discussion is limited to the relative range 
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between two orbiting bodies the issue of stability beeomes more eomplicated. The 
boundary region can typically be drawn such that Chase spacecraft stays within it. 
However, the open loop relative dynamics of spacecraft are inherently unstable, with 
non-zero initial conditions and external forces causing spacecraft to drift over time [70]. 
The A matrix of the linearized dynamics, shown in equation (4.16), yields six 
eigenvalues consisting of two repeated roots at zero and two repeated imaginary complex 
pairs. Therefore, the system is bounded-input, bounded-output (BIBO) unstable in the 
orbital plane and BIBO stable along the out-of-plane axis [71]. It is well known that 
close loop control can be used to achieve stable and convergent spacecraft relative 
motion. The stability of distributed spacecraft with individual feedback in a common 
reference field was rigorously addressed by Kang, Sparks, and Banda [72]. Although, 
their discussion was primarily based on formation maintenance it directly applies to 
convergence to any desired relative orbital position. The multiple spacecraft control 
algorithm balances the convergence to a desired goal position with the safety of collision 
avoidance. The convergence of the algorithm in a free environment can be explained by 
using common heuristic techniques and shown via simulation. However, in the presence 
of obstacles there may not be a unique solution. If the working space is not over 
constrained, multiple solutions will exist for the spacecraft path in the presence of 
multiple obstacles. There is also freedom in the closeness of the solution, based on 
acceptable goal achievement tolerance. Any attempt at a rigid stability proof may need to 
show convergence within a desired velocity and position ball, such as suggested by 
Neubauer and Swartwouf s research [55]. It is generally possible to select a closed map 
over which the algorithm is stable and robust. For these closure operations, preliminary 
initial conditions and simulation parameters must be carefully selected. 

Even if a global minimum exists, the spatial dynamics must allow freedom of 
motion for obstacle avoidance. Obstacles may be arranged, or maneuvered, to block the 
path of the spacecraft, beyond the capability of the control logic or the available control 
effort. It the absences of obstacles the multiple spacecraft control algorithm converges to 
the optimal trajectory. This convergence within a bounded area around an equilibrium 
state can be shown for each individual spacecraft, but can not be directly applied to the 
case of multiple spacecraft with uncertainty in relative states. Despite this lack of a 
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rigorous stability proof, the control algorithm can be shown to be rather robust in the 
presence of perturbations. The general stability of spacecraft formations in the presence 
of disturbances was shown by Acikmese et al [73]. The stability definition is motivated 
by BIBO definitions and the propagation of disturbances through the interconnected 
system. Sensor and actuator uncertainties can always be increased to a level that 
handicaps the control algorithm. Therefore, in the presence of high levels of noise and 
uncertainty, additional control logic needs to be implemented to ensure safety. The 
stability of control in the presence of multiple dynamic obstacles has yet to be proved. In 
this research global convergence was achieved by representing obstacles with spherical 
regions of influence. This minimizes the local minimization and chattering problems, 
since nonlinear dynamics act as persistent perturbations which will require control effort. 
Any local minimums encountered tend to be momentary and unstable. Also, only 
allowing obstacle influences to damp relative motion in the direction of the obstacle helps 
ensure stability. Motion along the path of the obstacle is restricted only to avoid 
collision, but movement along the convergent direction is unimpeded. This damping 
allows convergence as long as there is a free path to maneuver around an obstacle toward 
the goal. Any obstacle avoidance which allows divergent motion within the relative 
space must take into consideration the compounding effects of orbital rotations and 
escape velocities. 

Spacecraft translation actuation is often discretely implemented in the continuous 
dynamics of relative motion. It is worth mentioning that there is promising work in the 
area of hybrid control that lays the framework for evaluation of such systems [70] [74]. 
The hybrid control is well suited to the spacecraft translational problem with flow set 
(flow map) relating to free-floating continuous time and jump sets Oump rnaps) relating 
to control actuation at discrete time. The discussion involves weakly invariant sets and 
proper selection of Lyapunov functions. 

Both the stability and robustness of the close proximity control algorithm is 

demonstrated implicitly by Monte Carlo simulation results; refer to Chapter IX. There 

are three conditions which are considered control algorithm failures. These failures are 

based on the collision avoidance requirement, desire for short duration close proximity 

maneuvers, and the limited spacecraft propellant. First, any spacecraft collisions are 
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considered failures. Seeond, any maneuver durations lasting longer than 90 minutes is 
considered a failure. Finally, any spaeeeraft whieh uses all of its propellant is eonsidered 
a failure. None of these failures were detected during the Monte Carlo simulations. The 
elose proximity eontrol algorithm is considered stable due to the unfailing eonvergence to 
aceomplish the desired maneuver. Neither loeal minimums nor environmental 
disturbanees appeared to prevent Chase spaeeeraft from aehieving desired maneuver 
positions. Also, the elose proximity eontrol algorithm appears to be robust to a wide 
envelope of initial conditions and configurations. 
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VIII. CLOSE PROXIMITY OPERATIONS EVALUATION 


Control algorithm performance evaluation of a wide range of close proximity 
operations was conducted. For each maneuver, the multiple spacecraft are initially 
assumed to be within a spatial sphere with a one kilometer radius and negligible relative 
velocity to a Target spacecraft in a thrust free circular trajectory. For most close 
proximity maneuvers results cubic shaped spacecraft were simulated. The cubic shape 
proved to be the most challenging of the three basic shapes, spherical, cylindrical, and 
cubic spacecraft shapes for collision avoidance. More complicated shapes may require 
more conservative safety margins for the simultaneous maneuvering of multiple 
spacecraft. The communication of the goal, or target, position is the minimum 
information required. Desired mission ranges for inspection, regions of control for 
rendezvous, and spatial tolerances for docking may be dependent on both sensor and 
actuator performance. In particular collision avoidance is a balance between situational 
awareness, which is dependent upon sensor accuracy, and spacecraft’s reaction, which is 
based upon actuator capabilities. 

The primary parameters used to evaluate control algorithm performance are based 
on the maneuver duration time and required control effort. The duration of time required 
to successfully accomplish the desired close proximity maneuver is , in seconds. The 
general desire is that close proximity maneuvers can be conducted in approximately 30 
minutes. The control effort is related to the velocity change, Av in meters per second, 
required to complete the maneuver. Control effort should be efficient while maintaining 
desired performance. These two basic figures of merit, C and Av , must be accompanied 
by engineering analysis, since they are typically inversely related to each other. This is a 
conceptual simplification since these metrics are a result of the minimization of a cost or 
potential function with numerous constraints. Evaluation must take into account that any 
control effort which is severely, or continually, saturated may be a safety hazard in the 
vicinity of multiple spacecraft. Additional control effort may not be available to perform 
collision avoidance. Any control effort response that saturates the available control effort 
is denoted with an asterisk. Also, the overall convergence toward a goal should be 
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sufficiently damped in order to avoid oseillation in the goal region. Any oscillation in 
close vieinity to other spaeeeraft will result in pulsing responses of all other loeal 
spacecraft. 

Comprehensive performanee evaluation of the LQR/APF and APF eontrol 
algorithms was conducted for a wide range of close proximity operations, including 
convergence, rally, rendezvous, and doeking maneuvers. All simulated maneuvers 
ineluded orbital perturbations due to variations in the Earth's shape and mass, 
atmospherie drag on the spaeeeraft, third-body (Sun and Moon) forces, and solar- 
radiation pressure. Similarly, the gravity gradient, atmospheric drag, and solar pressure 
torques were included as perturbations in the attitude control loop. The close proximity 
operations evaluated all began with Chase spacecraft within 1.0 km of the goal position 
or Target spacecraft. The eonvergenee maneuver is simply moving to a goal position in 
free-space. The rally maneuver is gathering of multiple spacecraft to a common goal 
region in free spaee. Rendezvous maneuvers require the convergence of multiple 
spaeeeraft to a Target spaeeeraft. Doeking maneuvers require precise convergence to the 
outer boundary of a Target spacecraft while avoiding collision. The parameter results 
listed for eaeh maneuvers are representative of the general control algorithm 
performance. It should be understood that improvement to a particular parameter for a 
selected maneuver may be possible. The following results were generated using the same 
basie logic and gains for all maneuvers. 

A. CONVERGENCE MANEUVERS 

Close proximity operations begins when spacecraft independent eourse 
correetions and phasing has plaeed the spaeeeraft within a kilometer of eaeh other. The 
eonvergenee maneuver is intended to be a baseline maneuver, without collision 
avoidance, used for determining eontrol algorithm performance. In this maneuver, the 
Chase spaeeeraft maneuvers from its initial loeation to within 1.0 - 2.0 mm of goal 
position. Both the recursive LQR/APF and APF eontrol algorithms were successful in 
converging to within 1.0 mm of a goal position. The eonvergenee maneuver was used to 
tune algorithms gains for similar performanee based on maneuver duration and control 
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effort efficiency. This terminal range is much less than typically required for general 
close proximity control and serves to establish legitimacy for application of the control 
algorithm. 

The control algorithm performance results for six independent near and six 
independent far initial spacecraft positions are listed in Table 8.1 and Table 8.2, 
respectively. Each one of the near maneuvers starts approximately 100 meters way from 
a goal; whereas the far maneuvers start approximately 1,000 meters from a goal. 
Comparison with the two impulse maneuvers previously discussed, yields further insight 
into control algorithm performance. The first Chase spacecraft, listed in Table 8.1, 
arrives approximately 15 times faster with over 60 times more control effort than the two- 
impulse maneuver, which is discussed in Chapter IV.D. This is an extreme case of time 
and control effort trade-offs due to spacecraft phase dynamics. The second spacecraft, 
listed in Table 8.1, arrives approximately 3.7 times faster with only 1.5 times more 
control effort then the two-impulse maneuver, which is discussed in Chapter IV.D. The 
two-impulse maneuver is used for illustration purposes only, since execution of the two- 
impulse maneuver with precision necessary for close proximity operations is unfeasible. 
Realistic control response for multiple spacecraft close proximity maneuvers do not allow 
for implementation of two-impulse maneuvers. Comparison illustrates parameter trade¬ 
offs which must consider while developing multiple spacecraft control algorithms. 


Near Convergence Maneuver 

LQR/APF 

APF 

Near Convergence 

Initial RSW [0, 70, 0] m 

Av = 0.1877 m/s 

Av = 0.1905 m/s 

td= 1041 s 

td= 1264 s 

Near Convergence 

Initial RSW [50,-100, -50] m 

Av = 0.3105 m/s 

Av = 0.3123 m/s 

td= 1068 s 

td= 1298 s 

Near Convergence 

Initial RSW [100, 100, 100] m 

Av = 0.4900 m/s 

Av = 0.5121 m/s 

td= 1068 s 

td= 1317 s 

Near Convergence 

Initial RSW [100, 0, 0] m 

Av = 0.3077 m/s 

Av = 0.3215 m/s 

td= 1056 s 

td= 1284 s 

Near Convergence 

Initial RSW [-50, 100, -100] m 

Av = 0.3889 m/s 

Av = 0.3912 m/s 

td= 1082 s 

td= 1295 s 

Near Convergence 

Initial RSW [0, 0, 100] m 

Av = 0.2486 m/s 

Av = 0.2548 m/s 

td = 1062 s 

td= 1279 s 


Table 8.1 Six Spacecraft Near Convergence Maneuver. 
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For near convergence maneuvers the LQR/APF control algorithm is more control 
efficient and reaches the goal location faster then the APF control algorithm. This is as 
expected since the LQR/APF is not being affected by collision avoidance and is well 
within the linear performance region of the Clohessy-Wiltshire equations, refer to 
Equation (4.15). The APF has included a velocity ramping function to limit control 
saturation due to initial transients, refer to Equation (6.15). This allows for better 
comparison of performance in relative position, relative velocity, and control effort. The 
relative position, velocity and control effort performance of the APF control algorithm 
for the second Chase spacecraft, listed in Table 8.1, are shown in Figure 8.1 through 
Figure 8.3, respectively. The relative position, velocity and control effort performance of 
the EQR/APF control algorithm for the same Chase spacecraft are shown in Figure 8.4 
through Figure 8.6, respectively. 



Figure 8.1 Chase Spacecraft Relative Position Using APF for Near Convergence. 
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Figure 8.2 Chase Spaeeeraft Relative Veloeity Using APF for Near Convergenee. 



Figure 8.3 Chase Spacecraft Control Effort Using APF for Near Convergence. 
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Figure 8.4 Chase Spacecraft Relative Position Using LQR/APF for Near Convergence. 



Figure 8.5 Chase Spacecraft Relative Velocity Using LQR/APF for Near Convergence. 
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Figure 8.6 Chase Spacecraft Control Effort Using LQR/APF for Near Convergence. 

As shown, the LQR/APF and APF can both exhibit the desired close proximity 
performances. The control effort of the APF algorithm is primarily driven by the desired 
velocity. This is particularly apparent in the control effort response corresponding to 
velocity shaping function influence, as shown in Figure 8.3. By comparison, the 
LQR/APF control effort is balanced by the velocity and position states, as shown in 
Figure 8.6. The LQR/APF control algorithm generally results in an initial control action 
to begin convergence and a terminal control effort to stop at the desired goal. 

Additionally, far convergence maneuver results are listed in Table 8.2. The 
relative position, velocity and control effort performance of the APF control algorithm 
for the second Chase spacecraft, listed in Table 8.2, are shown in Figure 8.7 through 
Figure 8.9, respectively. The relative position, velocity and control effort performance of 
the LQR/APF control algorithm for the same Chase spacecraft are shown in Figure 8.10 
through Figure 8.12, respectively. The velocity shaping and control effort responses are 
more pronounced for these far maneuvers figures. Notice that higher velocity is achieved 
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during the execution of the far convergence maneuver. The APF control response 
reaches actuator saturation for the fourth Chase spacecraft, as denoted by the asterisk. 
This is generally an undesirable condition, but may not be completely avoidable in a 
dynamic environment without underutilizing the actuators. 


Far Convergence Maneuver 

LQR/APF 

APF 

Far Convergence 

Av = 2.5514 m/s 

Av = 2.2800 m/s 

Initial RSW [0, 1000, 0] m 

td= 1368 s 

td= 1446 s 

Far Convergence 

Av = 2.5454 m/s 

Av = 2.1884 m/s 

Initial RSW [412,-812,-412] m 

td= 1371 s 

td= 1453 s 

Far Convergence 

Av = 2.7375 m/s 

Av = 2.9295 m/s 

Initial RSW [575,575,575] m 

td= 1369 s 

td= 1459 s 

Far Convergence 

Av = 3.2204 m/s 

Av = 3.4209 m/s * 

Initial RSW [1000,0,0] m 

td= 1367 s 

td= 1449 s 

Far Convergence 

Av = 2.0063 m/s 

Av = 2.0201 m/s 

Initial RSW [0,0,1000] m 

td= 1389 s 

td= 1445 s 

Far Convergence 

Av = 2.9865 m/s 

Av = 3.231 m/s 

Initial RSW [707,707,0] m 

td= 1361 s 

td= 1454 s 

Table 8.2 Six Spacecraft Far Convergence M 

aneuver. 
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Figure 8.8 Chase Spacecraft Relative Velocity Using APF for Far Convergence. 
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Figure 8.10 Chase Spaeeeraft Relative Position Using LQR/APF for Far Convergenee. 
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Figure 8.11 Chase Spacecraft Relative Velocity Using LQR/APF for Far Convergence. 



Figure 8.12 Chase Spacecraft Control Effort Using LQR/APF for Far Convergence. 
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As expected, for both controllers the closer maneuvers took less time to complete. 
The LQR/APF results are generally more control efficient at similar, or faster, maneuver 
durations. For the closer rendezvous the multiple spacecraft LQR/APF was better in all 
metrics. This is primarily due to the LQR consideration of dynamics for the convergence 
control. In particular, the LQR algorithm results in a smooth initial increase in velocity 
while the APF requests a step increase in velocity. It was observed that the duration of 
some more distant maneuvers may be a few seconds shorter using the APF, if the APF 
initial velocity errors are not ramped. However, this is not considered significant, due to 
the average 30 minute maneuver duration. The 1.0 millimeter accuracy is used to show 
the capabilities of the algorithms, although sensor accuracy is not expected to be able to 
support this level of precision. Preliminary comparison to other APF applications, such 
as that presented in [52], [40] and [55], indicate that the algorithms accomplish higher 
precision within much shorter maneuver duration. Improved efficiency is generally 
accomplished, with direct comparisons strongly dependent on the density of the obstacle 
environment and initial conditions. Most previous algorithms delineate each phase of 
sequential maneuvering spacecraft with independent logic determining varying results, as 
in [54]. This makes the comparison of sequentially and simultaneous multiple spacecraft 
maneuvers in a dynamic environment difficult and unwarranted. However, this 
simultaneous and effective maneuvering of several spacecraft may be essential for a 
number of spacecraft missions. 

B. RALLY MANEUVERS 

The LQR/APF and APF control algorithms can be used to rally, or cluster, a 
group of spacecraft to a common point. Collision avoidance is used to avoid obstacles 
and other spacecraft while converging to within a desired range of a rally point. The 
results for three cubic Chase spacecraft simultaneously rallying proximity maneuver with 
collision avoidance are listed in Table 8.3. For these collision avoidance maneuvers, the 
goal position is a near rally point in free space. This requires the Chase spacecraft to 
converge to within 0.5 meters of the goal point while avoiding impact with other 
converging spacecraft. In addition, stationary obstacles are placed at positions along the 
path of each of the Chase spacecraft, for the sake of testing the algorithms. These 
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obstacles have an actual diameter of 2.0 meters and are positioned so that the Chase 
spacecraft encounter them at approximately 500-600 seconds into the maneuver. This 
point was selected to ensure that the obstaele was eneounter while the Chase spaceeraft 
was at relatively high veloeity. The three Chase spaeecraft motion during the rally 
maneuver is illustrated in Figure 8.13. Notice that the maneuver duration may lengthen 
for the later eongregating spacecraft, such as the case for the third Chase spacecraft listed 
in Table 8.3. This is due to other Chase spaeecraft arriving and loitering in the rally 
region. The spaceeraft whieh arrive first may limit the approaeh of the following 
spacecraft. This spatial constraint must be considered before determining the acceptable 
rally region in which the Chase spacecraft are being commanded. 



Table 8.3 Three Spaeecraft Near Rally Maneuver with Collision Avoidance. 
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Figure 8.13 Three Spaeeeraft Relative Motion During Rally Maneuver. 


The performance evaluation is once again extended to relatively far initial 
positions. The results for three cubic Chase spacecraft simultaneously rallying far 
maneuver with collision avoidance are listed in Table 8.4. 


Far Rally Maneuver 

LQR/APF 

APF 

Far Rally with Obstacle 

Av = 4.1887 m/s 

Av = 3.7436 m/s * 

Initial RSW [0, 1000, 0] m 

td= 1615 s 

td= 1243 s 

Far Rally with Obstacle 

Av = 4.0973 m/s * 

Av = 3.9324 m/s * 

Initial RSW [412, -812, -412] m 

td= 1503 s 

td= 1242 s 

Far Rally with Obstacle 

Av = 4.2541 m/s 

Av = 4.8076 m/s * 

Initial RSW [575, 575, 575] m 

td= 1487 s 

td= 1258 s 


Table 8.4 


Three Spacecraft Far Ra^ 


ly 


Maneuver with Collision Avoidance. 


In order to further challenge the control algorithms the simultaneous rally maneuvers 
were conducted with six Chase spacecraft. The near rally maneuver with collision 
avoidance results are listed in Table 8.5 and the far rally maneuver with collision 
avoidance results are listed in Table 8.6. In these six spacecraft maneuvers the rally 
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sphere is extended to 2.0 meters. Such each spacecraft must simultaneously converge to 
within 2.0 meters of the goal position while avoiding each other. The convergence 
sphere increases as the number of spacecraft increases, due to practical packing 
considerations. In realistic operational scenarios, safely margins would be maintained 
between spacecraft to allow for autonomous activity. 



Table 8.5 Six Spacecraft Near Rally Maneuver with Collision Avoidance. 



Table 8.6 Six Spacecraft Far Rally Maneuver with Collision Avoidance. 

The first Chase spacecraft to arrive in the goal region use less control effort. They 
do not need to avoid collisions with other spacecraft, whereas latter arriving spacecraft 
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may need to maneuver between to complete their rally maneuver. As the number of 
spacecraft increase direct comparison between spacecraft maneuvers become 
challenging. Each spacecraft maneuver is dependent upon every other spacecraft and 
obstacle within region of influence. Generally, the algorithms are comparable with 
maneuver durations being accomplished within a couple minutes of each other. 
Efficiency of the algorithms is dependent on the number and location of obstacles 
encountered and the desired relative velocities. Eor instance, the APE/EQR algorithm 
appears to be generally more efficient in a low density obstacle close proximity 
environment. However, the geometric based velocity management of the APE may 
converge faster in high density obstacle environments. This faster convergence comes at 
the risk of saturating the actuators. 

The EQR/APE and APE both continue to perform satisfactory in high density 
obstacle environments. However, the APE continues to saturate control effort. In some 
instances the saturation actually helps the APE control algorithms control efficiency, such 
as the first and second Chase spacecraft in Table 8.6. The efficiency appears to be better 
then the EQR/APE control algorithm, but this is only due to heavy actuator saturation 
which is generally undesired. The relative position, velocity and control effort 
performance of the APE control algorithm for the second Chase spacecraft, listed in 
Table 8.6, are shown in Eigure 8.14 through Eigure 8.16, respectively. The relative 
position, velocity and control effort performance of the EQR/APE control algorithm for 
the same Chase spacecraft are shown in Eigure 8.17 through Eigure 8.19, respectively. 
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Figure 8.15 Chase Spacecraft Relative Velocity Using APF for Far Rally. 
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Figure 8.16 Chase Spacecraft Control Effort Using APE for Ear Rally. 



Eigure 8.17 Chase Spacecraft Relative Position Using LQR/APE for Par Rally. 
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Figure 8.18 Chase Spacecraft Relative Velocity Using LQR/APF for Far Rally. 
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Figure 8.19 Chase Spacecraft Control Effort Using LQR/APF for Far Rally. 
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These results validate the LQR/APF eontrol algorithm for the multiple spaeecraft 
close proximity maneuvers. The saturation seen in the Chase spacecraft’s control effort 
due to the APF, refer to Figure 8.16, is not exhibited by the LQR/APF, refer to Figure 
8.19 The LQR/APF smoothly handles collision avoidance while successfully rallying 
multiple spacecraft to a desired location. 


C. RENDEZVOUS MANEUVERS 

Additionally, the LQR/APF and APF control algorithms with collision avoidance 
are used to avoid obstacles and other spacecraft while converging to within a specified 
range of a Target spacecraft’s outer boundary. For the three Chase spacecraft 
simultaneous rendezvous maneuvers the Chase spacecraft approach to within 5.0 
centimeters of the Target spacecrafts outer surface. The results for the three spacecraft 
simultaneously conducting a near rendezvous maneuver with collision avoidance are 
listed in Table 8.7. Similarly, the results for three spacecraft simultaneously conducting a 
far rendezvous maneuver with collision avoidance are listed in Table 8.8. For these 
collision avoidance maneuvers, the goal position is the center of the Target spacecraft. 
This requires that the Target spacecraft’s repulsion to allow the Chase spacecraft to 
converge while avoiding impact. Stationary obstacles are placed at positions along the 
unobstructed path of the Chase spacecraft. These obstacles are encountered at 
approximately 500-600 seconds into the maneuver, in order to ensure that the Chase 
spacecraft was at relatively high velocity. These obstacles have an actual diameter of 2.0 
meters. Generally, avoiding larger obstacles requires more control effort and time. The 
performance accuracy shows that convergence and rendezvous can be achieved while 
ensuring collision avoidance, as illustrated in Figure 8.20. Although, densely packed 
obstacle regions tend to result in a superposition of repulsion forces which keep later 
converging spacecraft at further distances away from the mutual goal location. 
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Near Rendezvous Maneuver 

LQR/APF 

APF 

Near with Obstacle 

Av = 0.3886 m/s 

Av = 0.3339 m/s 

Initial RSW [0, 70, 0] m 

td= 1342 s 

td= 1103 s 

Near with Obstacle 

Av = 0.6149 m/s 

Av = 0.5198 m/s 

Initial RSW [50,-100, -50] m 

td= 1314 s 

td= 1092 s 

Near with Obstacle 

Av = 0.8301 m/s 

Av = 0.8245 m/s 

Initial RSW [100, 100, 100] m 

td= 1540 s 

td=1117s 


Table 8.7 


Three Spaeeeraft Near Rendezvous Maneuver with Collision Avoidanee. 


Far Rendezvous Maneuver 

LQR/APF 

APF 

Far with Obstacle 

Initial RSW [0, 1000, 0] m 

Av = 4.1934 m/s 

Av = 3.7902 m/s * 

td= 1982 s 

td= 1276 s 

Far with Obstacle 

Initial RSW [412, -812, -412] m 

Av = 4.1455 m/s 

Av = 3.9995 m/s * 

td= 1720 s 

td= 1276 s 

Far with Obstacle 

Initial RSW [575, 575, 575] m 

Av = 4.4081 m/s 

Av = 4.8637 m/s * 

td= 1729 s 

td= 1287 s 


Table 8.8 Three Spacecraft Far Rendezvous Maneuver with Collision Avoidance. 



Figure 8.20 Three Spacecraft Collision Avoidance During Rendezvous Maneuver. 
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In order to further ehallenge the eontrol algorithms the simultaneous rendezvous 
maneuvers were eondueted with six Chase spaeeeraft. For the six Chase spaeeeraft 
simultaneous rendezvous maneuvers the Chase spaeeeraft approaeh to within 1.0 meter of 
the Target spaeeerafts outer surfaee. The near rendezvous maneuver with eollision 
avoidanee results are listed in Table 8.9 and the far rendezvous maneuver with collision 
avoidance results are listed in Table 8.10. 


Near Rendezvous Maneuver 

LQR/APF 

APF 

Near with Obstacle 

Initial RSW [0, 70, 0] m 

Av = 0.3138 m/s 

Av = 0.2763 m/s 

td= 1066 s 

td= 1011 s 

Near with Obstacle 

Initial RSW [50,-100, -50] m 

Av = 0.5294 m/s 

Av = 0.4424 m/s 

td= 1073 s 

td= 1017 s 

Near with Obstacle 

Initial RSW [100, 100, 100] m 

Av = 0.8097 m/s 

Av = 0.7362 m/s 

td= 1116s 

td= 1041 s 

Near with Obstacle 

Initial RSW [100, 0, 0] m 

Av = 0.4611 m/s 

Av = 0.5287 m/s 

td= 1049 s 

td= 1069 s 

Near with Obstacle 

Initial RSW [-50, 100, -100] m 

Av = 0.6939 m/s 

Av = 0.5854 m/s 

td= 1148 s 

td= 1026 s 

Near with Obstacle 

Initial RSW [0, 0, 100] m 

Av = 0.4605 m/s 

Av = 0.4551 m/s 

td= 1249 s 

td= 1159 s 


Table 8.9 Six Spacecraft Near Rendezvous Maneuver with Collision Avoidance. 


Far Rendezvous Maneuver 

LQR/APF 

APF 

Far with Obstacle 

Av = 4.1568 m/s 

Av = 3.7170 m/s * 

Initial RSW [0, 1000, 0] m 

td= 1504 s 

td= 1235 s 

Far with Obstacle 

Av = 4.0249 m/s 

Av = 3.9096 m/s * 

Initial RSW [412, -812, -412] m 

td= 1499 s 

td= 1225 s 

Far with Obstacle 

Av = 4.3766 m/s 

Av = 4.8191 m/s * 

Initial RSW [575, 575, 575] m 

td= 1523 s 

td= 1255 s 

Far with Obstacle 

Av = 4.8820 m/s * 

Av = 4.9273 m/s * 

Initial RSW [1000, 0, 0] m 

td= 1531 s 

td= 1251 s 

Far with Obstacle 

Av = 3.3261 m/s 

Av = 3.6335 m/s * 

Initial RSW [0, 0, 1000] m 

td= 1547 s 

td= 1232 s 

Far with Obstacle 

Av = 2.8049 m/s 

Av = 3.1474 m/s 

Initial RSW [707 ,707, 0] m 

td= 1317 s 

td= 1174 s 


Table 8.10 Six Spacecraft Far Rendezvous Maneuver with Collision Avoidance. 
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The efficiency of the maneuvers is not as clear in this maneuver due to non- 
optimal collision avoidance affecting both algorithms. The collision avoidance cause a 
decrease in acceleration and velocity as the Chase spacecraft approaches an obstacle. 
The LQR/APF algorithm tends to respond in a smooth manner while avoiding collisions. 
The APT algorithm tends to saturate the thrusters in the vicinity of obstacles. The rigid 
desired velocity following characteristic the APT algorithm may allow quick and efficient 
convergence during some maneuvers. However, there is a greater danger of collision as 
thrusters saturate and may not be able to respond to additional maneuver demands. The 
collision avoidance required at the Target spacecraft further decreases the efficiency of 
the LQR/APF control algorithm by making it stop short of its optimized goal. 

As expected, it generally takes longer for spacecraft to cover longer distance. 
Collision avoidance may further delay rendezvous due to the need to maneuver around 
obstacles. Notice that the duration of the simultaneous maneuvers increases for the latter 
arriving spacecraft due to repulsion of the other rendezvous spacecraft. The collision 
avoidance algorithm logic ensures that converging spacecraft are not perturbed by latter 
converging spacecraft. This ensures safety in the convergence, but may cause the 
rendezvous of later spacecraft to be delayed due to the congestion at the shared goal 
position. For instance, the Chase spacecraft avoiding obstacles can be delayed due to the 
obstacle avoidance maneuver in free space and experience additional convergence delay 
due to latter arrival at the rendezvous. This last delay usually affects the third, or more, 
spacecraft approaching the rendezvous point. This delay can be resolved by dedicating 
different rendezvous points for each spacecraft. This concept is explicitly be addressed 
by assigning different goal locations for each spacecraft, such as in docking maneuvers. 

D. DOCKING MANEUVERS 

The final stage and ultimate goal of rendezvous may be the docking of multiple 

spacecraft. The two spacecraft docking maneuver is the basis for on orbit servicing and 

assembly. As multiple spacecraft are required to perform docking maneuvers, several 

potential complications arise. The docking mechanisms and the docking order need to be 

addressed. The forces and torque tolerance of the docking mechanism and the overall 

spacecraft need to be considered. Also, the docking mechanisms must be arranged on 
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each spacecraft to allow for clear fields of view for sensors and clear approach zones. 
For spacecraft the docking order may be predetermined. This is typically the case for 
heterogeneous spacecraft that must be assembled in a specific order or manner. For 
homogenous spacecraft, the order of docking may not be as important, but there may be 
limitations due to docking mechanism number and position. Each docking mechanism 
may be function as a female, male, or both. Also, the spacecraft structure may also limit 
the docking mechanism placement on the spacecraft. For instance, a possible cubic 
spacecraft may dock on any of its six sides; refer to Figure 8.21. The dotted lines 
represent male to female docking orientations. A spacecraft with only one male and 
female connection is very limited in versatility of assembly scenarios. 



Figure 8.21 Cubic Spacecraft Docking Positions. 

The developed LQR/APF and APF control algorithms were evaluated for the 
docking maneuver. As with the other maneuvers, the docking maneuver time and fuel 
efficiency are the primary metrics used for evaluating the performance of a control 
algorithm. Each of Chase spacecraft was assigned a desired goal location on the outer 
surface of the Target spacecraft. Both EQR/APE and APE control algorithms with 

collision avoidance are used to avoid obstacles and other spacecraft while converging to 
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within 2.0 millimeter of the doeking position on the Target spaceeraft’s outer boundary. 
The results of the three spacecraft near docking proximity maneuver with collision 
avoidance are listed in Table 8.11. The same initial positions and stationary obstacles 
were selected for direct comparison. The results of the three spacecraft far docking 
proximity maneuver with collision avoidance are listed in Table 8.12. 


Near Docking Maneuver 

LQR/APF 

APF 

Near with Obstacle 

Initial RSW [0, 70, 0] m 

Av = 0.3355 m/s 

Av = 0.2959 m/s 

td= 1135 s 

td= 1282 s 

Near with Obstacle 

Initial RSW [50,-100, -50] m 

Av = 0.5694 m/s 

Av = 0.5397 m/s 

td= 1118s 

td= 1296 s 

Near with Obstacle 

Initial RSW [100, 100, 100] m 

Av = 0.9915 m/s 

Av = 0.9823 m/s * 

td= 1355 s 

td= 1386 s 


Table 8.11 Three Spacecraft Near Docking Maneuver with Collision Avoidance. 


Tab 


Far Docking Maneuver 

LQR/APF 

APF 

Far with Obstacle 

Av = 4.1575 m/s 

Av = 3.7561 m/s * 

Initial RSW [0, 1000, 0] m 

td= 1522 s 

td= 1477 s 

Far with Obstacle 

Av = 4.4056 m/s 

Av = 4.0500 m/s * 

Initial RSW [412, -812, -412] m 

td=2119s 

td= 1479 s 

Far with Obstacle 

Av = 4.6149 m/s 

Av = 5.003 m/s * 

Initial RSW [575, 575, 575] m 

td= 1878 s 

td= 1531 s 


e 8.12 Three Spacecraft Far Docking Maneuver with Collision Avoidance. 


For sub-centimeter docking precision, the exact relative velocity of the docking 
port must be taken into consideration. Objects appearing to be stationary in relative 
position actually have a velocity which directly relates to the orbital rotation of the RSW 
frame. This velocity offset can limit precession if not determined for all objects 
appearing relatively stationary in the RSW frame. This velocity offset, , is 

determined by simply multiplying the docking port’s offset position from the Target 
center by the orbital rotation rate. 
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This resulting velocity is that which is required for an object to appear stationary 
with respect to the center of the RSW frame. Without this correction for LEO orbits, the 
resulting convergence error is approximately 3.0-4.0 cm for every meter of offset. This 
offset velocity correction is only necessary for relatively fixed objects and does not 
impact freely orbiting objects. 

Ultimately the docking mechanism must be sized for the spacecraft. 
Consideration must be made for forces and torques due to manner of docking, such as 
grappling, magnetic, or male/female. The final position and velocity error control 
algorithm must be within the tolerance of the docking mechanism. For instance, if the 
final position is expected to be within 2 cm then the docking mechanism may need to 
have an adaptor with a radius of 2 cm plus the size of the docking lever; refer to Figure 
8.22. The velocity error of the final position will result in forces and torques imparted 
between spacecraft as they make contact. This docking impact may result in transient 
translational and rotations which need to be controlled. For instance, even if a mother 
spacecraft can support the docking of multiple spacecraft, the docking order may need to 
be staggered in order to allow the assembled space structure to stabilize after each 
docking. This is less significant for small spacecraft docking to relatively more massive 
spacecraft. 
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The control algorithms are applied to the more complicated maneuver involving 
six Chase spacecraft simultaneously docking to a single Target spacecraft. For the six 
Chase spacecraft docking maneuvers, the ports are centered on each side of a cubic 
Target spacecraft at RSW locations of [1, 0, 0], [0, -1, 0], [0, 1, 0], [-1, 0, 0], [0, 1, 0], and 
[0, -1, 0], respectively. The docking locations were randomly assigned and not related to 
the initial positions of the Chase spacecraft. This allows for robustness evaluation of the 
control algorithm, since it requires that each Chase spacecraft maneuver in close 
proximity. Both LQR/APF and APF control algorithms with collision avoidance are used 
to avoid obstacles and other spacecraft while converging to within 2.0 millimeter of their 
assigned docking position on the Target spacecraft’s outer boundary. The results for six 
Chase spacecraft conducting simultaneously near docking maneuvers, with collision 
avoidance, are listed in Table 8.13. The results for six Chase spacecraft simultaneously 
conducting far docking maneuvers, with collision avoidance, are listed in Table 8.14. 


Near Docking Maneuver 

LQR/APF 

APF 

Near with Obstacle 

Av = 0.3355 m/s 

Av = 0.2959 m/s 

Initial RSW [0, 70, 0] m 

td= 1135 s 

td= 1282 s 

Near with Obstacle 

Av = 0.5694 m/s 

Av = 0.5648 m/s * 

Initial RSW [50,-100, -50] m 

td= 1118 s 

td= 1304 s 

Near with Obstacle 

Av = 0.9915 m/s 

Av = 1.1742 m/s * 

Initial RSW [100, 100, 100] m 

td= 1355 s 

td= 1516 s 

Near with Obstacle 

Av = 0.6902 m/s 

Av = 0.7264 m/s * 

Initial RSW [100, 0, 0] m 

td= 1307 s 

td= 1444 s 

Near with Obstacle 

Av = 0.9969 m/s 

Av = 0.7238 m/s * 

Initial RSW [-50, 100, -100] m 

td= 1521 s 

td= 1397 s 

Near with Obstacle 

Av = 0.7399 m/s 

Av = 0.6944 m/s * 

Initial RSW [0, 0, 100] m 

td= 1794 s 

td= 1923 s 


Table 8.13 Six Spacecraft Near Docking Maneuver with Collision Avoidance. 
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Far Docking Maneuver 

LQR/APF 

APF 

Far with Obstacle 

Av = 4.1792 m/s 

Av = 3.7513 m/s * 

Initial RSW [0, 1000, 0] m 

td= 1530 s 

td= 1478 s 

Far with Obstacle 

Av = 4.5243 m/s 

Av = 4.1084 m/s * 

Initial RSW [412, -812, -412] m 

td= 1870 s 

td= 1488 s 

Far with Obstacle 

Av = 4.7083 m/s 

Av = 5.2832 m/s * 

Initial RSW [575,575, 575] m 

td= 1719 s 

td= 1549 s 

Far with Obstacle 

Av =4.9268 m/s * 

Av = 5.2024 m/s * 

Initial RSW [1000, 0, 0] m 

td= 1520 s 

td= 1602 s 

Far with Obstacle 

Av = 3.6151 m/s 

Av = 3.8136 m/s * 

Initial RSW [0, 0, 1000] m 

td= 1678 s 

td= 1496 s 

Far with Obstacle 

Av = 3.0789 m/s 

Av = 5.1804 m/s * 

Initial RSW [707, 707, 0] m 

td= 1463 s 

td= 1509 s 


Table 8.14 


Six Spacecraft Far Docking Maneuver with Co 


lision Avoidance. 


The relative position, velocity and control effort performance of the APF control 
algorithm for the second Chase spacecraft, listed in Table 8.14, are shown in Figure 8.23 
through Figure 8.25, respectively. The relative position, velocity and control effort 
performance of the LQR/APF control algorithm for the same Chase spacecraft are shown 
in Figure 8.26 through Figure 8.28, respectively. These responses are similar to those 
seen for previous maneuvers, which illustrate the general similarity in multiple spacecraft 
close proximity maneuvers. Although, additional control response is often needed in the 
later stages of precision docking maneuvers, refer to the control effort response in Figure 
8.28. 
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Figure 8.23 Chase Spaceeraft Relative Position Using APF for Far Rally. 



Figure 8.24 Chase Spacecraft Relative Velocity Using APF for Docking. 
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Figure 8.25 Chase Spaeeeraft Control Effort Using APE for Ear Booking. 



Figure 8.26 Chase Spaeeeraft Relative Position Using LQR/APF for Far Docking. 
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Figure 8.27 Chase Spaeeeraft Relative Veloeity Using LQR/APF for Far Docking. 



Figure 8.28 Chase Spacecraft Control Effort Using LQR/APF for Far Docking. 
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The most interesting part of multiple spaeeeraft docking is best shown with 
animated graphic. This allows evaluation of control logic and convergence performance 
with respect to how each spacecraft maneuvers with relationship to both time and space. 
Without animation, a basic assessment of the LQR/APF and APF controller’s relative 
position performance can be illustrated with a path comparison. The relative path of the 
second Chase spacecraft, from Table 8.13 and Table 8.14, is shown in Figure 8.29. The 
commanded path for each of the control algorithms is obviously different, with variations 
having cascading effects for collision avoidance with other spacecraft. The stationary 
obstacle collision avoidance and docking regions are highlighted in Figure 8.30 and 
Figure 8.31, respectively. The stationary obstacle avoidance response is straight forward; 
however the docking regions path requires further explanation. The Chase spacecraft 
arrives at the Target spacecraft in a different order based on LQR/APF or APF control. 
Using the APF, this Chase spacecraft arrives as the second spacecraft in the docking 
region, so modest collision avoidance in the docking region is required. For the 
LQR/APF this particular Chase spacecraft arrives last and must fully employ collision 
avoidance. Variations in commanded path, arriving sequence, and relative position of the 
other Chase spacecraft make one to one performance comparison challenging. The 
control behavior of each Chase spacecraft may appear to be intelligent behavior, but is 
actually the result of the control algorithm’s basic computation and logic. These control 
algorithms can be expanded to handle additional logical situations of particular 
maneuvers. 
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R-axis, m 



Figure 8.29 LQR/APF and APF Path Comparison. 



S-axis, m 

Figure 8.30 LQR/APF and APF Collision Avoidanee Path Comparison. 
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Figure 8.31 LQR/APF and APF Docking Region Path Comparison. 

Selecting which spacecraft should dock at which docking mechanism may be a 
function of a higher order control algorithm. For instance a docking assignment 
algorithm may select the closest relative spacecraft to each docking mechanism that is 
available. Multiple spacecraft docking assignment depends on the mission and 
characteristic of each spacecraft. One example of a candidate task assignment algorithm 
was developed in [1] and could be adapted for the docking task. This algorithm is based 
on dynamics programming and allows for the spacecraft order to be based on fuel 
efficiency of the maneuver. This algorithm can be readily adapted to potential function 
relative positioning without a centralized computation and a simple communication 
protocol. However, once a hierarchical tasking algorithm is applied additional 
communication and computation may be required. 
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E. MANEUVER EVALUATION CONCLUSIONS 

Based on the wide range of elose proximity maneuvers evaluated, the LQR/APF 
eontrol algorithm is an exeellent candidate for multiple spacecraft close proximity 
operations. The convergence rate and control effort efficiency is comparable to that of a 
highly tuned APF control algorithms. The incorporation of the linearized dynamics in the 
LQR/APF control algorithm allows it to overcome limitations of APF control algorithms 
for actual implementation. The iterative calculations allow for practical computation for 
real time applications. The LQR/APF simulations demonstrate successful completion of 
multiple spacecraft conducting simultaneous, rally, rendezvous, and docking maneuvers. 
Inclusion of high fidelity spacecraft perturbation forces during the simulations proves that 
the LQR/APF algorithm is also effective in disturbance rejection. 

The control effort expended during the close proximity maneuver can be related 
to the quantity of propellant used. Spacecraft have an estimated margin of 3-6% for 
propellant mass, refer to Table 3.2. Based on 100 kg spacecraft, the initial propellant 
mass would be at least 3.0 kg. The far simultaneous docking maneuvers are potential the 
most control effort intensive maneuvers. The control effort for the six Chase spacecraft 
docking maneuvers, refer to Table 8.14, are represented in propellant mass in Figure 
8.32. From this perspective, a Chase spacecraft could perform about ten far docking 
maneuvers. This is an extremely conservative estimate based on far initial positions, 
short maneuver duration, modest volume of propellant, and a dense obstacle 
environment. The average Chase should be able to perform several close proximity 
maneuvers. Therefore, the LQR/APF algorithm appears to be a promising new 
development for the field of multiple spacecraft close proximity maneuver control. 
Monte-Carlo method analysis allows for reasonable statistical estimates of the mean and 
standard deviation of maneuver C and Av . Refer to Chapter IX for multiple spacecraft 
control algorithm Monte Carlos simulation results. 
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Spacecraft Mass, kg 



Figure 8.32 Propellant Usage in Relation to Mass. 
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IX. MONTE-CARLO ANALYSIS OL CLOSE PROXIMITY 
MULTIPLE SPACECRALT MANEUVERS 


Numerous simulations of multiple spaeecraft close proximity maneuvers were 
conducted in order to generate a sample distribution of maneuver parameters. Using 
Monte-Carlo methods, estimates of the mean and standard deviation of maneuver 
duration, t^, and delta-v, Av, were determined. Two hundred convergence, rally, 

rendezvous, and docking maneuver simulations were conducted for both the APF and 
LQR/APF control algorithms. Each simulation involved six Chase spacecraft performing 
simultaneous maneuvers. The Law of large numbers permits the approximation of 
sample statistics via Monte-Carlo methods. As expected, the large sample size generally 
approaches a Gaussian distribution. The normalized data distribution allows for 
estimates of the maneuver parameter means and standard deviations. The statistical data 
is presented in both per spacecraft and per maneuver format. The per spacecraft statistics 
use each Chase spacecraft of each maneuver for a total sample size of 1,200 spacecraft. 
The per maneuver statistics use the maximum parameters of each maneuver for a total 
sample size of 200. The average close proximity maneuver duration is valuable to 
spacecraft operators. Similarly, the average Av is valuable to both spacecraft designers 
and mission planners. 

None of the three control algorithm failures conditions, as discussed in Chapter 
VIED, were experienced by either the refined APF or the developed LQR/APF control 
algorithm. First, no spacecraft collisions were detected. Second, all spacecraft 
maneuvers were successfully performed within 90 minutes. Finally, no spacecraft was 
required to use all of its propellant during maneuvering. Therefore, both control 
algorithms proved to be effective in performing close proximity operations. A statistical 
analysis of both the APF and LQR/APF during the convergence and docking maneuvers 
follows. 
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A. INITIAL CONDITIONS AND SPACECRAFT PARAMETERS 

For each simulation, the Chase spacecraft initial positions were uniformly 
randomly distributed while all other simulation parameters were maintained. Each of the 
six Chase spacecraft was initially positioned within a 1.0 km sphere with respect to a 
Target spacecraft position. However, the Chase spacecraft initial position was assumed 
to be at least 10 m from the Target. This stand off range, representing the center-to- 
center distance between the Target and Chase spacecraft, was used so that spacecraft size 
variation of a few meters could be readily simulated. The Target spacecraft was assumed 
to be in a circular LEO of 500 km altitude. The initial relative velocity was assumed to 
be negligible. This neutral velocity state in the relative frame, suggests an elliptical 
orbital phase for the Chase spacecraft. Eor the sake of control evaluation, this neutral 
situation is reasonable and serves to avoid bias due to favorable velocity conditions. The 
relatively high velocity management of both the APE and LQR/APE allows for some 
initial velocity. The initial velocity was used as an experimental control variable, while 
the initial position is treated as a uniformly randomly distributed independent variable. 

The initial range of each Chase spacecraft is within a 10 - 1,000 m sphere of the 
Target Spacecraft. The initial positions of all Chase spacecraft, with respect to the Target 
spacecraft, are shown in Eigure 9.1. This includes all 1,200 Chase spacecraft of the 200 
simulations. The initial position of the first Chase spacecraft of each simulation is shown 
in Eigure 9.2. This subset of 200 spacecraft initial positions shows that the Chase 
spacecraft are randomly distributed for each simulation. The Chase spacecraft Monte- 
Carlo simulation initial range statistics are listed in Table 9.1. Eirst, the initial range 
mean and standard deviation is listed for all 1,200 Chase spacecraft in the 200 
simulations. The Chase spacecraft are uniformly randomly distributed, as shown in 
Eigure 9.3. Next, the mean and standard deviation is listed for the maximum initial 
Chase spacecraft range of each of the 200 simulations. The maximum Chase spacecraft’s 
initial range distribution for each maneuver, as shown in Eigure 9.4. This maximum 
initial range metric drives the overall and Av for each multiple spacecraft maneuver. 

Both the LQR/APE and APE control algorithm are analyzed over the same random range 
distribution. This allows for direct comparison of performance for each maneuver. 
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Figure 9.1 Initial Position of All Chase Spaeeeraft. 



Figure 9.2 Initial Position of First Chase Spacecraft. 


Initial Range Statistics 

Mean 

Standard Deviation 

Initial Range (1200 samples) 

501.8 m 

282.1 m 

Max Initial Range (200 samples) 

840.0 m 

122.0 m 


Table 9.1 Chase Spacecraft Range Distribution Statistics. 
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Figure 9.3 Initial Chase Spacecraft Range Distribution. 
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Maneuver Max Initial Range, m 

Figure 9.4 Maximum Maneuver Initial Range Distribution. 
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The Chase spaeecraft parameters were held eonstant for all simulations. All 
spaceeraft are homogeneous cubes of 1.0 m. Each spacecraft has an initial mass of 100 
kg; however the mass of each varies as it uses propellant. The general Monte-Carlo 
simulation method parameters are listed in Table 9.2. Each high fidelity spacecraft 
model is subject to all of the perturbing forces and torques discussed in Chapter III. 


Monte Carlo Analysis Simulation Parameters 

Number of Simulations 

200 

Spacecraft Shape 

Cubic 

Spacecraft Size 

1.0 m 

Initial Spacecraft Mass 

100 kg 

Initial Propellant Mass 

3 kg 

Target Spacecraft Orbit 

Circular orbit at 500 km altitude 

Number of Chase Spacecraft per Maneuver 

6 

Chase Spacecraft Initial Position 

10 - 1,000 m from Target 

Chase Spacecraft Initial Velocity 

Negligible relative velocity 

Spacecraft Perturbation Eorces 

High Eidelity Model 


Table 9.2 Monte Carlo Simulation Parameters. 


B. MONTE-CARLO ANALYSIS OF CONVERGENCE MANEUVERS 

The simple convergence maneuver, once again, serves as the baseline for 
evaluating control algorithm performance. The individual spacecraft distribution 
statistics for C and Av are listed in Table 9.3. Overall, the EQR/APE performs simple 

convergence maneuvers better, on a per spacecraft basis. The maneuvers are 
accomplished in less time using less control effort. This is as expected due to the 
iterative optimization of the EQR algorithm. The EQR/APE’s and APE’s convergence Q 
spacecraft distributions are shown in Eigure 9.5 and Eigure 9.6. The EQR/APE’s and 
APE’s convergence Av spacecraft distributions are shown in Eigure 9.7 and Eigure 9.8. 
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Convergence Spacecraft Statistics 
(1,200 Samples) 

LQR/APF 

APF 

td 

Mean 

1212.9 s 

1364.0 s 

Standard Deviation 

102.6 s 

56.5 s 

Av 

Mean 

1.252 m/s 

1.276 m/s * 

Standard Deviation 

0.708 m/s 

0.732 m/s 


Table 9.3 Convergence Spacecraft Statistics. 
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Figure 9.5 LQR/APF Spacecraft Convergence Maneuver Duration Distribution. 
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Figure 9.6 APF Spacecraft Convergence Maneuver Duration Distribution. 
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Figure 9.7 LQR/APF Spacecraft Convergence Delta-V Distribution. 
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Figure 9.8 APF Spacecraft Convergence Delta-V Distribution. 


The statistical analysis is extended to a per maneuver basis. The maximum C and 

Av are selected for each six Chase spacecraft of a convergence maneuver. Additionally, 
the Av of all six maneuvering Chase spacecraft are summed to obtain the total Av for 
each maneuver. The mean and standard deviation of these metrics are listed in Table 9.4. 
The LQR/APF control algorithm performs simple convergence maneuvers well. The 
maneuver durations are shorter on average. Both the maximum and total maneuver Av 
have a smaller average and tighter standard deviation. The tighter standard deviation is 
helpful for spacecraft designers and mission planners. The LQR/APF’s and APF’s 
convergence maneuver maximum C distributions are shown in Figure 9.9 and Figure 

9.10. The maximum duration for each convergence maneuver is similar to the maximum 
initial range distribution. The LQR/APF’s and APF’s convergence maneuver maximum 
Av distributions are shown in Figure 9.11 and Figure 9.12. The LQR/APF’s and APF’s 
total convergence maneuver Av distributions are shown in Figure 9.13 and Figure 9.14. 
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Convergenc 

(2( 

e Maneuver Statistics 
)0 Samples) 

LQR/APF 

APF 

Max td 

Mean 

1330.7 s 

1423.4 s 

Standard Deviation 

38.7 s 

17.5 s 

Max Av 

Mean 

2.147 m/s 

2.219 m/s 

Standard Deviation 

0.404 m/s 

0.459 m/s 

Total Av 

Mean 

7.515 m/s 

7.656 m/s 

Standard Deviation 

1.899 m/s 

1.951 m/s 


Table 9.4 Convergence Maneuver Statistics. 
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Figure 9.9 LQR/APF Convergence Maneuver Maximum Duration Distribution. 
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Figure 9.10 APF Convergence Maneuver Maximum Duration Distribution. 



Figure 9.11 LQR/APF Convergence Maneuver Maximum Delta-V Distribution. 
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Figure 9.12 APF Convergence Maneuver Maximum Delta-V Distribution. 



Figure 9.13 LQR/APF Convergence Maneuver Total Delta-V Distribution. 


147 
































































22 -1-1-!-1-r 



Total Maneuver Delta-V, m/s 

Figure 9.14 APF Convergence Maneuver Total Delta-V Distribution. 

During convergence maneuvers, the APF control algorithm significantly saturated 
the control actuation in approximately 3-6% of all convergence maneuvers. These would 
represent situations where no further control effort is available for collision avoidance. 
This saturation allows the APF control effort to appear comparatively efficient, but is not 
a desired state for high density obstacle environments. 


C. MONTE-CARLO ANALYSIS OF RALLY MANEUVERS 

For the Monte-Carlo six spacecraft rally maneuvers the rally sphere is 3.5 meters. 
Each cubic spacecraft must simultaneously converge to within 3.5 meters of the goal 
position while avoiding each other. This spherical range represents the maximum cross- 
section of two cubic spacecraft. The individual spacecraft distribution statistics for C 
and Av are listed in Table 9.5. Overall, the LQR/APF performs six spacecraft rally 
maneuvers better, on a per spacecraft basis. The maneuvers are accomplished in slightly 
more time using less control effort. This is as expected due to the strict velocity control 
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of the APF control algorithm. The LQR/APF’s and APF’s rally spacecraft 

distributions are shown in Figure 9.15 and Figure 9.16. The LQR/APF has a couple 
duration outliers due to Chase spacecraft loitering in front of later arriving Chase 
spacecraft. The LQR/APF’s and APF’s rally Av spacecraft distributions are shown in 
Figure 9.17 and Figure 9.18. 


Rally Spacecraft Statistics 
(1,200 Samples) 

LQR/APF 

APF 

td 

Mean 

1170.5 s 

1065.2 s 

Standard Deviation 

179.1 s 

90.2 s 

Av 

Mean 

1.151 m/s 

1.167 m/s 

Standard Deviation 

0.706 m/s 

0.720 m/s 


Table 9.5 Rally Spacecraft Statistics. 



Spacecraft Maneuver Duration, s 


Figure 9.15 LQR/APF Spacecraft Rally Maneuver Duration Distribution. 
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APF Spacecraft Rally Maneuver Duration Distribution. 
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Figure 9.17 LQR/APF Spacecraft Rally Delta-V Distribution. 
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Spacecraft Delta-V, m/s 

Figure 9.18 APF Spacecraft Rally Delta-V Distribution. 


The statistical analysis is extended to a per rally maneuver basis, with the mean 
and standard deviation of and Av metrics listed in Table 9.6. The LQR/APF control 
algorithm performs rally maneuvers well. The average maneuver duration is slightly 
longer, but both the maximum and total maneuver Av has smaller average and tighter 
standard deviation. The APF control algorithm significantly saturated the control 
actuation in approximately 5-10% of all rally maneuvers. The LQR/APF’s and APF’s 
rally maneuver maximum distributions are shown in Figure 9.19 and Figure 9.20. The 

maximum duration for each rally maneuver is similar to the maximum initial range 
distribution. The LQR/APF’s and APF’s rally maneuver maximum Av distributions are 
shown in Figure 9.21 and Figure 9.22. The LQR/APF’s and APF’s total rally 
maneuver Av distributions are shown in Figure 9.23 and Figure 9.24. 
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Rally M 

(2( 

aneuver Statistics 
)0 Samples) 

LQR/APF 

APF 

Max td 

Mean 

1341.9 s 

1143.4 s 

Standard Deviation 

184.1 s 

21.2 s 

Max Av 

Mean 

2.055 m/s 

2.109 m/s 

Standard Deviation 

0.389 m/s 

0.454 m/s 

Total Av 

Mean 

6.907 m/s 

7.003 m/s 

Standard Deviation 

1.808 m/s 

1.881 m/s 


Table 9.6 Rally Maneuver Statistics. 



Maximum Maneuver Duration, s 


Figure 9.19 LQR/APF Rally Maneuver Maximum Duration Distribution. 
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9.20 APF Rally Maneuver Maximum Duration Distribution. 
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Figure 9.21 LQR/APF Rally Maneuver Maximum Delta-V Distribution. 
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Figure 9.22 APF Rally Maneuver Maximum Delta-V Distribution. 
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Figure 9.23 LQR/APF Rally Maneuver Total Delta-V Distribution. 
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Figure 9.24 APF Rally Maneuver Total Delta-V Distribution. 


D. MONTE-CARLO ANALYSIS OF RENDEZVOUS MANEUVERS 

The six spaceeraft rendezvous maneuver sphere is similar to that of the rally 
maneuver. Eaeh spaceeraft must simultaneously converge to within 3 meters of the 
surface of the Target spacecraft while avoiding each other. This allows rendezvous 
within twice the maximum cross-section of the cubic spacecraft. Any closer rendezvous 
requires thorough operational consideration of practical maneuver space and end goal 
points. The individual spacecraft distribution statistics for and Av are listed in Table 

9.7. Similar to the rally maneuver, the LQR/APF performs six spacecraft rendezvous 
maneuvers in slightly more time using less control effort. The LQR/APF’s and APF’s 
rendezvous Q spacecraft distributions are shown in Figure 9.25 and Figure 9.26. As seen 
in the rally maneuver, the LQR/APF has a couple duration outliers. The LQR/APF’s and 
APF’s rendezvous Av spacecraft distributions are shown in Figure 9.27 and Figure 9.28. 
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Rendezvous Spacecraft Statistics 
(1,200 Samples) 

LQR/APF 

APF 

td 

Mean 

1171.9s 

1065.0 s 

Standard Deviation 

181.3 s 

90.1 s 

Av 

Mean 

1.150 m/s 

1.166 m/s 

Standard Deviation 

0.704 m/s 

0.719 m/s 


Table 9.7 Rendezvous Spaeeeraft Statisties. 



Spacecraft Maneuver Duration, s 


Figure 9.25 LQR/APF Spaeeeraft Rendezvous Maneuver Duration Distribution. 
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Figure 9.26 APF Spacecraft Rendezvous Maneuver Duration Distribution. 
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Figure 9.27 LQR/APF Spacecraft Rendezvous Delta-V Distribution. 
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Figure 9.28 APF Spacecraft Rendezvous Delta-V Distribution. 


The statistical analysis is extended to a per rendezvous maneuver basis, with the 
means and standard deviations of and Av listed in Table 9.8. The LQR/APF control 
algorithm performance is similar to that of the rally maneuver, with slightly longer 
duration and smaller average Av and tighter Av standard deviation. Again, the APF 
control algorithm significantly saturated the control actuation in approximately 5-10% of 
all rendezvous maneuvers. The LQR/APF’s and APF’s rendezvous maneuver maximum 
Q distributions are shown in Figure 9.29 and Figure 9.30. The maximum duration for 

each rendezvous maneuver is similar to the maximum initial range distribution. The 
LQR/APF’s and APF’s rendezvous maneuver maximum Av distributions are shown in 
Figure 9.31 and Figure 9.32. The LQR/APF’s and APF’s total rendezvous maneuver Av 
distributions are shown in Figure 9.33 and Figure 9.34. 
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Rendezvous 

(2( 

Maneuver Statistics 
)0 Samples) 

LQR/APF 

APF 

Max td 

Mean 

1347.8 s 

1143.2 s 

Standard Deviation 

188.0 s 

21.3 s 

Max Av 

Mean 

2.054 m/s 

2.107 m/s 

Standard Deviation 

0.386 m/s 

0.454 m/s 

Total Av 

Mean 

6.902 m/s 

6.994 m/s 

Standard Deviation 

1.793 m/s 

1.878 m/s 


Table 9.8 Rendezvous Maneuver Statistics. 
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Figure 9.29 LQR/APF Rendezvous Maneuver Maximum Duration Distribution. 


159 





































25 -!-1-1-1-r 



1080 1100 1120 1140 1160 1180 1200 

Maximum Maneuver Duration, s 

Figure 9.30 APF Rendezvous Maneuver Maximum Duration Distribution. 
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Figure 9.31 LQR/APF Rendezvous Maneuver Maximum Delta-V Distribution. 
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Figure 9.32 APF Rendezvous Maneuver Maximum Delta-V Distribution. 
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Figure 9.33 LQR/APF Rendezvous Maneuver Total Delta-V Distribution. 
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Figure 9.34 APF Rendezvous Maneuver Total Delta-V Distribution. 

E. MONTE-CARLO ANALYSIS OF DOCKING MANEUVERS 

The six spacecraft simultaneous docking maneuver requires that all six spacecraft 
avoid each other while converging to within 2.0 millimeter of their assigned docking 
position on the Target spacecraft’s outer boundary. The docking port positions, on the 
center of each cubic face of the Target spacecraft, are randomly assigned. The individual 
spacecraft distribution statistics for and Av are listed in Table 9.9. The LQR/APF still 

performs well, however the collision avoidance maneuvering close to the Target favors 
the APF controller. The maneuvers are accomplished in slightly more time and more 
control effort. This is as expected due to the intrinsic collision avoidance capability of 
the APF control algorithm. Although, the APF control efficiency is only achieved by 
saturating the available control actuation. The LQR/APF’s and APF’s docking 
spacecraft distributions are shown in Figure 9.35 and Figure 9.36. The LQR/APF’s and 
APF’s docking Av spacecraft distributions are shown in Figure 9.37 and Figure 9.38. 
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Docking Spacecraft Statistics 
(1,200 Samples) 

LQR/APF 

APF 

td 

Mean 

1460.9 s 

1438.0 s 

Standard Deviation 

202.6 s 

122.0 s 

Av 

Mean 

1.423 m/s 

1.382 m/s 

Standard Deviation 

0.735 m/s 

0.764 m/s 


Table 9.9 Docking Spacecraft Statistics. 



Spacecraft Maneuver Duration, s 


Figure 9.35 LQR/APF Spacecraft Docking Maneuver Duration Distribution. 
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Figure 9.36 APF Spacecraft Docking Maneuver Duration Distribution. 
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Figure 9.37 LQR/APF Spacecraft Docking Delta-V Distribution. 
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Figure 9.38 APF Spacecraft Docking Delta-V Distribution. 


The statistical analysis is extended to a per docking maneuver basis, with the 
means and standard deviations of and Av listed in Table 9.10 The LQR/APF control 
algorithm performance continues to maintain a tighter Av standard deviation. This 
confirms the LQR/APF control algorithm’s more predictable performance. While on the 
other hand, the APF control algorithm significantly saturated the control actuation in 
approximately 10-20% of all docking maneuvers. This increase in actuator saturation is 
undesirable and increases the risk of collision in high density obstacle regions. The 
LQR/APF’s and APF’s docking maneuver maximum Q distributions are shown in Figure 
9.39 and Figure 9.40. The LQR/APF’s and APF’s docking maneuver maximum Av 
distributions are shown in Figure 9.41 and Figure 9.42. The LQR/APF’s maximum Av 
is lower on average with a tighter distribution. The LQR/APF’s and APF’s total docking 
maneuver Av distributions are shown in Figure 9.43 and Figure 9.44. 
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Docking ]> 
(2( 

klaneuver Statistics 
)0 Samples) 

LQR/APF 

APF 

Max td 

Mean 

1672.0 s 

1602.2 s 

Standard Deviation 

194.5 s 

162.9 s 

Max Av 

Mean 

2.348 m/s 

2.362 m/s 

Standard Deviation 

0.394 m/s 

0.456 m/s 

Total Av 

Mean 

8.537 m/s 

8.291 m/s 

Standard Deviation 

1.918 m/s 

1.989 m/s 


Table 9.10 Docking Maneuver Statistics. 
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Figure 9.39 LQR/APF Docking Maneuver Maximum Duration Distribution. 
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Figure 9.40 APF Docking Maneuver Maximum Duration Distribution. 
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Figure 9.41 LQR/APF Docking Maneuver Maximum Delta-V Distribution. 
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Figure 9.42 APF Docking Maneuver Maximum Delta-V Distribution. 
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Figure 9.43 LQR/APF Docking Maneuver Total Delta-V Distribution. 
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Figure 9.44 APF Docking Maneuver Total Delta-V Distribution. 


F. MONTE-CARLO ANALYSIS CONCLUSIONS 

The Monte-Carlo method analysis of close proximity maneuvers confirms that the 
LQR/APF control algorithm is a practical candidate for multiple spacecraft close 
proximity operations. Estimates of the mean and standard deviation of maneuver and 

Av show that average proximity maneuver control effort, Av, efficiency is generally 
better than that of a highly tuned APF control algorithm. The LQR/APF showed a 0.5- 
1.0% efficiency improvement on a per spacecraft maneuver basis. The standard 
deviation of the LQR/APF control effort is consistently narrower then that of the APF. 
The LQR/APF showed a 10-20% narrower standard deviation. This narrow Av standard 
deviation is valuable to both spacecraft designers and mission planners. It allows 
effective propellant sizing for close proximity operations. It also gives operational 
planners a useful tool for developing and forecasting maneuvers. The average mission 
duration, Q, of all close proximity operations were maintained below the desired 30 

minutes. The wider Q standard deviation of the LQR/APF control algorithm is due to 
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velocity variations allowed by the algorithm. This standard deviation of approximately 
three minutes is acceptable for spacecraft operations which have traditionally been 
measured is terms of hours or days. Both the LQR/APF’s and APF’s Av efficiency 
improves if the convergence rate of the Chase spacecraft is slowed. Therefore, the 
relative velocity can be used as a design trade-off between maneuver duration and 
efficiency. 

Based on this Monte-Carlo analysis, the LQR/APF control algorithm appears 
suitable for application to emerging multiple spacecraft operations. Both the control 
efficiency and maneuver duration are reasonable for current spacecraft designs. Based on 
this analysis, the average Chase spacecraft could perform 20 - 40 close proximity 
maneuvers. This is more maneuvering than typically discussed in even the most 
aggressive spacecraft servicing operations. Therefore, the LQR/APF algorithm appears 
to be practical for multiple spacecraft close proximity maneuver control. Variation in 
spacecraft physical characteristics and orbital assumptions may cause some fluctuations 
in the total number of close proximity maneuvers which can be performed. However, it 
the LQR/APF control algorithm performs reliably for a wide range of maneuvers. 
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X. PERFORMANCE VALIDATION AND VISUALIZATION 


A refined and validated multiple spaceeraft six DOF simulation was configured 
for the purpose of control algorithm development during close proximity operations. The 
simulator incorporates a six DOF MATLAB and Simulink numerical model with 3D 
STK visualization. Both the model validation and 3D visualization are crucial to 
successful engineering evaluation during control algorithm development. A MATLAB- 
STK simulation interface was developed in order to validate the spacecraft models and 
provide detailed animation of simulations [75]. 

An effective test scenario is one which dependably simulates the environment in 
which the control algorithm is expected to operate. The application of the control 
algorithm for use on multiple spacecraft in close proximity operations drives the 
requirements that it be tested with computer-generated orbital dynamics and kinematics. 
Spacecraft model validation gives confidence that results are consistent and reliable. 
This research produced a method of spacecraft model validation by comparison with STK 
spacecraft analysis software developed by AGI [36]. STK is used as an orbital 
propagator for both simulation and emulation of the desired spacecraft models. 

Accurate rendering of 3D spacecraft during compound close proximity maneuvers 
permit additional confidence in derived results. In addition to spacecraft model 
validation, STK can be used for detailed animation of spacecraft simulations. 2D 
animations based on complex numerical translational and attitude data often fail to 
convey true physical relationships. Even with multiple animated views these 2D 
representations may lead to conceptual misunderstandings. By utilizing 3D simulations 
more accurate engineering analysis can be conducted, allowing for undesirable 
performance to be discovered and corrected. This visualization is especially important 
for missions which require simultaneous control of multiple spacecraft maneuvering in 
close proximity. STK visualization can enhance spacecraft engineering analysis of 
relative translational motion and attitude while allowing for variations in spacecraft 
parameters and constraints. An overview of the MATLAB-STK simulation interface 
used for model validation and visualization is presented. 
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A. MATLAB-STK SIMULATION INTERFACE 

Due to its ubiquitousness in engineering applieations, a MATLAB based interfaee 
was desired. The common understanding of MATLAB code among both electrical and 
mechanical engineering researchers allows for a common starting point for multiple 
spacecraft simulation. In order to develop high fidelity spacecraft models, MATLAB and 
Simulink code is extremely useful and relevant. The high level language allows for ease 
in defining multiple parameters and numerical computation. The STK core modules 
allow for analytical and numerical orbital propagation of multiple spacecraft. Additional 
educational modules allow STK to be used in conjunction with MATLAB [39]. For this 
development MATLAB version 7.3.0267 (R2006b) and STK version 7.0.1 were utilized 
[39] [36]. Instructions for the installation and configuration of STK and MATLAB can 
be found at the AGI website [36]. By transferring information between MATLAB and 
STK, MATLAB constructed simulation dynamics and kinematics can be verified and 
visualized. 

The MATLAB-STK interface takes advantage of STK’s standard connection 
protocol for initializing STK from MATLAB. STK developers have enabled their code 
to be executed by MATLAB via an application program interface (API) called 
AgConnect [36]. The connection allows for properly formatted commands to be sent via 
a Transmission Control Protocol and Internet Protocol (TCP/IP) port. By exploiting this 
feature, formatted data, such as ephemeris and attitude date can be passed between 
applications. For instance, the interface allows all spacecraft physical characteristics and 
simulation parameters to be defined in MATLAB and related to STK. Also, STK data 
can be passed to MATLAB for analysis The MATLAB-STK simulation interface 
exploits the strength of each software application. The general functions and information 
flow between MATLAB and STK are shown in Fig. 1. This is a broad overview of the 
primary functions and how both applications are employed. 
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Figure 10.1 MATLAB-STK Simulation Interface Overview. 

1. Overview of MATLAB/Simulink Spacecraft Modeling 

Generalized spacecraft characteristics are modeled using MATLAB and Simulink 
[39]. The MATLAB and Simulink software components were modular structured 
allowing for variation in the number of spacecraft, orbital perturbations, spacecraft and 
obstacles parameters, and desired maneuvers. A high level Simulink spacecraft model is 
shown in Figure 10.2. The primary blocks are labeled with their functional calculations. 
The labeled blocks may have several sub layers with corresponding MATLAB code. The 
wiring has been simplified for visual flow, and some variables are shared between the 
kinematics and dynamics blocks. All multiple spacecraft simulation parameters are 
assigned, or defined, in MATLAB and implemented by a single Simuli nk model. 
Spacecraft states are concatenated into vectors and passed through the Simulink model. 
This vector of multiple spacecraft states passing through a single Simulink model is more 
computationally efficient and flexible then having multiple Simulink models to run 
simultaneously. The modeled space environment perturbations included variations in the 
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Earth's shape and mass, atmospheric drag on the spacecraft, third-body (Sun and Moon) 
forces, and solar-radiation pressure, and mass variation due to thruster firings; refer to 
Chapter III.E. 
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Eigure 10.2 High Eevel Simulink Model for Multiple Spacecraft Simulation. 


Several parameters must be defined, or assigned, before multiple spacecraft 
simulations can be executed. These parameters can include the number of simulations, 
initial time, selection of perturbations, number of spacecraft and obstacles, selection of 
control algorithms, and duration of simulation. Eor each simulation, a reference start time 
must be selected. This time can either be the current clock time, or a default reference 
time, which will be used for perturbation calculations, such as solar drag and third body 
effects. In this research, the Elniversal Time (UTl) in the format of [year, month, day, 
hour, minute, second] was used. With respect to this time, it is typical to select a duration 
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and time step for the simulation. A subset of simulation conditions are summarized in 
Table 10.1. The defining of spacecraft proximity operation to be within one kilometer is 
consistent with other researchers [76]. Although, this research considers rendezvous a 
subset of the close proximity maneuvers within the one kilometer range. It is worth 
mentioning that this meaning may vary from the phased rendezvous and proximity 
operations mission sequences outlined in the STK/Astrogator module [36]. 


Target Spacecraft 

Target Minimum Altitude 

300 km 

Target Maximum Altitude 

2,000 km 

Chase Spacecraft 

Number of Chase Spacecraft 

1 - 14 

Chase Spacecraft 
Initial Position 

R-axis 

10- 1,000 m 

S-axis 

10- 1,000 m 

W-axis 

10- 1000 m 

Chase Spacecraft 
Initial Velocity 

R-axis 

0 m/s 

S-axis 

0 m/s 

W-axis 

0 m/s 


Table 10.1 Close Proximity Spacecraft Maneuver Simulation Parameters. 


In addition to simulation parameters, there are adjustable parameters for each 
spacecraft, such as initial position, attitude, size, shape, mass, and control actuators. 
Refer to Chapter III.G. for discussion of these parameters. The initial position of 
spacecraft could be specifically assigned, as in Chapter VIII, or randomly determined, as 
in Chapter IX. Selection of the desired maneuver prompts the fine-tuning of the control 
algorithm logic, such as decreasing the region of influence due to other Chase spacecraft 
during the terminal stage of docking. Visualization of the system responses required 3D 
evaluations of multiple spacecraft operating in close proximity. Although, 3D 
visualization is possible using MATLAB it is limited and leads to some programming 
challenges. For instance, changing views during animations may be beyond the ability of 
most users. For a majority of applications simple 3D point mass dynamic representation 
with 2D plots of associated metrics might be acceptable, such as used in [58]. However, 
close proximity maneuvers of high fidelity 3D spacecraft require more realistic 


175 



simulation for detailed performance evaluation of collision avoidance and docking. The 
visualization of spacecraft orbits and maneuvers is a specific field which is greatly 
simplified by the use of STK. 

2. Overview of Satellite Tool Kit (STK) 

The STK is a useful spacecraft environment simulator and propagator with a 
sophisticated graphical user interfaces (GUI) for spacecraft particular applications. 
However, since STK is not as commonly used in the engineering fields as MATLAB, 
development of a simulation interface between the two proved be useful. In order to 
visualize and evaluate the controller’s performance, the MATLAB generated data was 
passed to STK. The spacecraft modeling and visualization capabilities of STK and its 
related modules and plug-ins are extensive. The visual analysis of the multiple spacecraft 
is invaluable during the evaluation of complex maneuvers. Both attitude and 
translational dynamics can be viewed from a wide variety of coordinate systems. Sensor 
field-of-view can be added to spacecraft models to provide an additional level of fidelity 
to the simulations. The STK/Connect module provides the means for STK to 
communicate with applications, such as MATLAB, through the use of a TCP/IP socket. 
This allows data to be transferred between MATLAB and STK. The sending and 
receiving of information to and from STK is intended to accurately model sensor 
performance between multiple spacecraft. Additional development may use the 
STK/Communications module to emulate communication propagation delay, frequency 
and bandwidth limitation, and bit error rates between all spacecraft [77]. This 
communication evaluation could be used to verily that higher order controller commands 
are properly implemented. 

Typical STK users select scenario parameters via the GUI interface. This allows 
primary objects, such as spacecraft, to be loaded with desired constraints. Additionally 
sensors can be defined and attached to the satellites. The GUI has numerous drop-down 
menus with many layers for defining basic characteristic, 2D graphic animation, 3D 
graphic animation, and related constraints. A sample STK scenario GUI interface is 
shown in Figure 10.3. The numerous button toolbars are listed along the top of the GUI. 
The primary objects and attachments are listed along the left column. The 3D graphics 
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animation is shown in the scenario screen, with the general STK user page and 2D 
graphic animation tabbed along the bottom of the screen. 



Figure 10.3 Sample STK Scenario GUI. 

For those interested, several additional STK supported modules may be useful for 
spacecraft proximity operations research. These modules supplement the core STK 
visualization environment and assist in the evaluation of mission sequences. In 
particular, the Astrogator module can be used to support detailed maneuver analysis and 
operations [76]. The standard API links directly to the STK simulation environment. 
The Astrogator module is particularly useful for detailed mission analysis bring more 
fidelity to the entire mission control phasing and sequences. Further integration and 
discussion concerning the Astrogator module is beyond the scope of this research, refer to 
[76] and [36]. 
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Similarly, the Advanced Close Approach Tool (ACAT) may be useful in collision 
avoidance maneuver planning. This STK module allows for additional situational 
awareness of any spacecraft or object in a referenced database. The ACAT module 
enables proximity indications and visual cues for variable close approach calculations. 
The use of uncertainty ellipsoids is more of a collision prediction tool for spacecraft path 
planning; refer to [47]. The analysis may be useful when evaluating potential sensor 
failure and communication link limitations. These modules may be more useful as space 
mission analysis and situational awareness applications. 

3. MATLAB-STK Interface 

The MATLAB-STK interface allows overall simulation parameters and spacecraft 
physical characteristics to be defined in MATLAB and related to STK via mexConnect 
commands, or formatted native STK commands. This interface technique takes advantage 
of STK’s standard connection protocol for initializing STK from MATLAB. Once 
initialized, native STK commands can be called from MATLAB. The specific 
requirements of each particular STK command will determine which path and parameters 
are necessary. In addition to commands, formatted ephemeris and attitude files can be 
passed from MATLAB to STK. Likewise, data can be pulled from STK by using the 
stkReport command. Pulling ephemeris from STK and passing it to the MATLAB 
engine allows for comparison of independently generated STK and MATLAB spacecraft 
propagation. Evaluating the results of this comparison allows for validation of a 
developed MATLAB model based on the STK’s High Precision Orbital Propagator 
(HPOP). The satellite of interest must be initialized and assigned with HPOP, before the 
propagation parameters and perturbations can be tailored via HPOP commands. By 
systematically enabling and comparing perturbation effects, each component of the 
spacecraft model can be validated. Once the numerical model has been validated, 
complex multiple spacecraft maneuvering can be animated for further assessment. 
However, precision close proximity maneuvers require refinement of the typical STK 
animation and views. For instance, the spacecraft’s physical characteristics, such as size, 
shape and mass, may need to be modified. These parameters can be assigned as 
MATLAB variables and executed in STK by placing them as parameters in STK native 
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commands. The STK graphical model is primarily modified by using an assortment of 
STK Visual Option (VO) commands. The MATLAB-STK simulation interface code is 
written in MATLAB {.m file format) with some MATLAB/Simulink files nested within 
it. These sub-files serve as modular components of the MATLAB-STK interface, which 
allow for easy simulation variation and modification between users. 


a. Instructions to Establish MATLAB-STK Interface Configuration 
For the standard MATLAB-STK connection, both MATLAB and STK 
applications should be loaded and launched. Once launched, STK can be initialized from 
MATLAB by using the command stkinit or agilnit. Information on the path 
setting established can be found by using agiGetConfig. The next step in the 
connection processes is to open a socket, by using stkOpen. This configures the path 
and assigns a socket variable and a MATLAB variable, called STKError, for reference. 
Via the established MATLAB/STK connection, STK can be commanded by using either 
a select number of mexConnect commands or native STK commands. The command 
paths for these two methods are slightly different. The mexConnect commands are a 
limited subset of core MATLAB/STK interface commands which can be found by 
exploring MATLAB Help. These mexConnect commands are all prefixed with stk, 
such as the stkinit commands mentioned above. Since the mexConnect commands 
are limited, there is a general command which allows native STK commands to be 
executed from MATLAB. The general execution of STK commands via MATLAB can 
be conducted by using stkExec (ConID, 'Command Path Parameter' ). A 
full listing of native STK commands which can be implemented in this fashion are listed 
in the STK Help menu under the path <Automate/Extend/Integrate>, <Command 
Listings>, <Alphabetical Listing>. The first step is usually to open a new STK scenario 
and ensure that any previously scenarios are closed. This ensures that a clean STK 
workspace is being established. The initial MATLAB-STK interface code is as follows: 

stkinit 

%initializes the STK/MATLAB Interface 
remMachine = stkDefaultHost; 

delete(get(0 children ')) ; %clear open MATLAB charts 

conID=stkOpen(remMachine); %Open the Connect to STK 
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%%Check to see if a scenario is open 
scen_open=stkValidScen; 
if scen_open == 1 

rtn = questdlg (' Close the current STK Scenario?'); 
if -strcmp(rtn, 'Yes' ) 
stkClose(conID) 
return 
else 

stkUnloadV ('/*') 

end 

end 

scen_nam=[ 'Multi_Spacecraft' ,num2str(CONST.simnum)]; 
stkNewObj( '/' , 'Scenario' ,scen_nam); 


The con ID variable serves as a numeric representation of the established connection. 
This conID variable will be used often in stkExec commands. The STK scenario 
name was determined by a MATLAB variable, such as CONST. simnum, which can 
selected by the user to represent the number of simulations or scenarios desired. The 
sceneries name was then used to establish the new STK scenario. The stkNewObj 
command path is establishing the new scenario at the highest level path level. All 
subsequent STK objects, such as spacecraft, will be assigned under this current scenario. 


b. STK Scenario Time Synchronization 

The selected simulation initial date and time must be properly formatted 
and passed to STK. The date and time can be represented in several formats. One such 
format is [year, month, day, hour, minute, second], with a specific example to [2007, 07, 
21, 12, 30, 0]. However, for STK applications the date and time must be rearranged and 
passed in the format of [21 Jul 2007 12:30:00.0]. If the MATLAB datestr command 
is used to determine the date and time, then a method of re-formatting and passing the 
date and time to STK is as follows: 


para_now=datestr(now); 
para_now(3)= ' ' ; 

para_now(7)= ' ' ; 

PARA.TIME=clock; %sets formatted time for MATLAB 

stkSetTimePeriod(para_now,para_now, 'GREGUTC') ; 
stkSetEpoch(para_now, 'GREGUTC' ); 
stkSyncEpoch; 

newpara=strcat (' SetValues_ para_now, o.2 0.1'); 

newpara(10)= ' ' ; 

rtn=stkConnect(conID, 'Animate' ,[ 'Scenario/' ,scen_nam], newpara); 
rtn=stkConnect (conID, 'Animate' , [ 'Scenario/' ,scen_nam], 'Reset' ); 
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In the newpara variable, the seeond to last number is animation time step, in seconds, 
and the last number is the highest speed or refresh rate, in seconds. The scale of these 
rates will influence the simulation by changing the animation of the scenario. These 
values can also be assigned as variables. When variables are passed into the body of an 
STK formatted command, brackets must be employed. This is apparent in the 
stkConnect command, where the variable scen_nam is used. Code listed in this 
research is only intended to be representative of the capability of the MATLAB-STK 
simulation interface, since there may be numerous methods of producing similar results. 


c. Simple Satellite Object 

Multiple spacecraft simulations require several spacecraft to be created. 
These spacecraft will need to be assigned physical and mass characteristics, attitude, and 
initial position and velocity. Along with the time determined above, these will serve as 
the basic parameters for initializing orbital propagation. The following sample code 
shows a spacecraft, called Target, being created; 


stkNewObj ('*/', 'Satellite ' , 'Target' ); 

stkExec(conID, 'VO */Satellite/Target Model Filename"s/c model"'); 
stkExec (conID, ['SetMass */Satellite/Target Value',... 
num2str(Mass)]) 

StkExec (conID, ['SetMass */Satellite/Target Matrix',... 
num2str(Inertia)]); 

StkExec(conID, [ 'SetAttitude */Satellite/Target Profile InertFix 
Quat' , num2str(QuaternionO), ' "CentralBody/Earth J2000"']); 
stkSetPropCart ( ' */Satellite/Target' , ' HPOP ' , ' J2000 ' , tstart, ... 
tstopdummy,stepsize,orbitepoch,STATE.ri(1:3)',STATE.vi(1:3)'); 


The VO command is used to call a general spacecraft graphical model which is available 

to STK. The spacecraft graphical model can be selected from a default list, usually 

located at C : \Program Files\AGI\STK 7\STKData\VO\Models\Space, or 

a custom user generated spacecraft graphical model. The development of simple 

spacecraft graphical models will be discussed in more details in Chapter X.C. In this 

sample code several variables are passed from MATLAB to STK. First, the mass 

characteristics are described by the scalar Mass and the vector Inertia. Inertia is 

the inertia matrix of the spacecraft in vector form. Next, the initial quaternion attitude of 

the Target spacecraft, called QuaternionO, is assigned. For this research, an inertial 

fixed attitude about the Earth is assigned using quaternion in the inertial fixed coordinate 
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system (J2000). The propagation requires an initial start and stop time, called tstart 
and tstopdummy respectively. The stop time is temporarily assigned at this point to 
allow for assignment of the HPOP integrator. Next, the orb it epoch variable is 
assigned. This variable is usually equal to zero. Finally, the initial position vector, called 
STATE, ri, and velocity vector, called STATE, vi, are passed into STK. These 
parameters define the initial Target spacecraft states, but additional HPOP settings will 
allow for selection of integration and perturbations features. These will enable for our 
systematic spacecraft numerical model validation. 

B. SPACECRAFT MODEL VERIFICATION 

The high fidelity six DOF spacecraft model and dynamics was verified by 
comparison of developed MATLAB and Simulink orbital modeling with custom STK 
propagators. The STK High Precision Orbital Propagator (HPOP) was used to ensure 
spacecraft propagation conditions matched all model variations. Due to its commonality 
in both Simulink and STK, fourth order Runga-Kutta numerical integration method was 
used. Details on numerical integration applied to orbital propagation are provided in 
[32]. The position and velocity states of non-maneuvering spacecraft were compared for 
different propagation variations. The modeled space environment perturbations, 
including variations in the Earth's shape and mass (J2-J4 coefficients), atmospheric drag 
on the spacecraft, third-body (Sun and Moon) forces, and solar-radiation pressure, were 
sequentially compared. In addition, the simple two-body dynamics serves as a baseline 
case and all of these perturbations were included into a total perturbation case. The 
comparisons of the MATLAB/Simulink and STK propagations show that spacecraft 
model validation via a MATLAB-STK simulation interface is achievable. 

I. MATLAB-STK Model Validation Interface 

The multiple spacecraft model was modularly configured to allow each of the 
perturbation models to be sequentially evaluated and validated. The modular 
perturbations are triggered by an enabling constant. For each module a variable, called 
CONST. PERT. choice, enabled the proper Simuli nk blocks and assigned the 
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corresponding HPOP settings for STK. The MATLAB and Simulink numerical 
modeling of the perturbations was conducted as discussed in Chapter III and [32], The 
desired precision of the numerical representations can vary as preferred by the user. 


a. HPOP Integrator Setting 

The HPOP command allows for assignment of both the numerical 
integration method, such as fourth order Runge-Kutta, and the Earth Gravitational model, 
such as EGM96 [36][32]. The HPOP command is of the general form, 
stkExec (ConID, 'HPOP Path Parameter' ) . Once the HPOP is assign to a 
spacecraft the perturbation forces can be further defined. The desired perturbations can 
be related to the enabling condition on the Simulink model, so that the various 
perturbations can be independently evaluated. The HPOP settings can be initialized for 
each spacecraft with the following code; 

StkExec(conID, 'HPOP */Satellite/Target Integrator IntegMethod RK4 
ReportOnFixedStep On'); 

if CONST.PERT.choice==l; %Earth Oblateness 

StkExec(conID, 'HPOP */Satellite/Target Force Gravity "C:/Program 
Files/AGI/STK 7/STKData/CentralBodies/Earth/EGM96.grv" 4 O'); 
StkExec(conID, 'HPOP */Satellite/Target Drag Off'); 

StkExec(conID, 'HPOP */Satellite/Target Force SolarRad Off'); 

StkExec(conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Moon 
Off' ); 

StkExec(conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Sun 
Off' ); 

elseif CONST.PERT.choice==2 ; %Aero Drag 

StkExec(conID, 'HPOP */Satellite/Target Force Gravity "C:/Program 
Files/AGI/STK 7/STKData/CentralBodies/Earth/EGM96.grv" 0 O'); 
StkExec(conID, ['HPOP */Satellite/Target Drag On ' ,Cd, ' 0.01 "1976 
Standard"' ]); 

StkExec(conID, 'HPOP */Satellite/Target Force SolarRad Off'); 

StkExec (conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Moon 
Off' ); 

StkExec(conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Sun 
Off' ); 

elseif CONST.PERT.choice==3; %Solar Drag 

StkExec(conID, 'HPOP */Satellite/Target Integrator IntegMethod RK4 
ReportOnFixedStep On'); 

StkExec(conID, 'HPOP */Satellite/Target Force Gravity "C:/Program 
Files/AGI/STK 7/STKData/CentralBodies/Earth/EGM96.grv" 0 O'); 
StkExec(conID, 'HPOP */Satellite/Target Drag Off'); 

StkExec(conID, ['HPOP */Satellite/Target Force SolarRad On ' ,Cd, ' 
0.01'] ) ; 

StkExec(conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Moon 
Off' ); 
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stkExec(conID, 'HPOP */ 
Off' ); 

elseif CONST.PERT.choice= 
StkExec(conID, 'HPOP */ 
ReportOnFixedStep On') 
StkExec(conID, 'HPOP */ 
Files/AGI/STK 7/STKDat 
StkExec(conID, 'HPOP */ 
StkExec(conID, 'HPOP */ 
StkExec(conID, 

On FromCB' ); 

StkExec(conID, 

On FromCB' ); 

elseif CONST.PERT.choice= 
StkExec(conID, 'HPOP */ 
ReportOnFixedStep On') 
StkExec(conID, 'HPOP */ 
Files/AGI/STK 7/STKDat 
StkExec(conID, ['HPOP * 
Standard"' ]); 

StkExec(conID, ['HPOP * 
0 . 01 '] ); 

StkExec(conID, 'HPOP */ 
On FromCB' ); 

StkExec(conID, 

On FromCB' ); 
else 

StkExec(conID, 'HPOP */ 
ReportOnFixedStep On') 
StkExec(conID, 'HPOP */ 
Files/AGI/STK 7/STKDat 
StkExec(conID, 


Satellite/Target Force ThirdBodyGravity Sun 
==4; %3-Body 

Satellite/Target Integrator IntegMethod RK4 

r 

Satellite/Target Force Gravity "C:/Program 
a/CentralBodies/Earth/EGM96.grv" 0 O'); 
Satellite/Target Drag Off'); 

Satellite/Target Force SolarRad Off'); 

HPOP */Satellite/Target Force ThirdBodyGravity Moon 

HPOP */Satellite/Target Force ThirdBodyGravity Sun 

==5; %Total 

Satellite/Target Integrator IntegMethod RK4 

r 

Satellite/Target Force Gravity "C:/Program 
a/CentralBodies/Earth/EGM96.grv" 4 O'); 
/Satellite/Target Drag On ' ,Cd, ' 0.01 "1976 

/Satellite/Target Force SolarRad On ' ,Cd, ' 

Satellite/Target Force ThirdBodyGravity Moon 

HPOP */Satellite/Target Force ThirdBodyGravity Sun 

%Two Body with no perturbations 
Satellite/Target Integrator IntegMethod RK4 


Satellite/Target Force Gravity "C:/Program 
a/CentralBodies/Earth/EGM96.grv" 0 O'); 

HPOP */Satellite/Target Drag Off'); 

StkExec(conID, 'HPOP */Satellite/Target Force SolarRad Off'); 

HPOP */Satellite/Target Force ThirdBodyGravity Moon 


StkExec(conID, 

Off' ); 

StkExec(conID, 'HPOP */Satellite/Target Force ThirdBodyGravity Sun 
Off' ); 
end 


Selection of the perturbation parameter, CONST . PERT . choice, activates the desired 
perturbation code. The selection may include all or none of the considered perturbations. 
For instance, the selection of CONST . PERT . choice=6 results in the simple two-body 
propagation of the spacecraft orbit. Additionally, the two numbers in the Gravity 
settings command, at the beginning of each conditional setting, define the degree and 
order of the zonal coefficient terms, respectively [36]. The coefficients can be 
synchronized with any particular user defined model by modifying the .grv file selected. 
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As previously discussed, variable can be passed into the STK command. For instance, 
the coefficient of drag, represented by Cd, was incorporated into the Drag and 
SolarRad command settings. 

Once the HPOP settings are completed the actual propagation can be 
computed. The initial stop time, tstopdummy, is replaced with the final desired 
propagation time, t stop. The spacecraft propagation command is as follows: 

stkSetPropCart ( ' */Satellite/Target' , ' HPOP ' , ' J2000 tst art, t stop,... 
steps1ze,orbitepoch,STATE.ri(1:3) ',STATE.vi(1:3) ') ; 

where the propagation parameters, stepsize, orbitepoch, STATE, ri, and 
STATE . vi are the same as in the initial propagation command; refer to Chapter X.A.3.C. 


b. Output STK Data to MATLAB Workspace 
The data from STK propagation can be passed to the MATLAB 
Workspace. This will allow any desired post processing and evaluation of the data. 
Proper variable assignment of the data allows for state variable comparison. STK 
spacecraft state data can be passed in several formats. For this research the ECI position 
and velocity were desired. Sample code for passing spacecraft ephemeris is as follows: 


[stkData,stkName]=stkReport( '*/Satellite/Target' , 'J2000 ECI 
Position Velocity'); 

STATE.STK.time=stkFindData(stkData{1}, 'Time' ); 

STATE.STK.r(:,1)=stkFindData(stkData{l}, 'x' ) ; 

STATE.STK.r(:,2)=stkFindData(stkData{1}, 'y ' ) ; 

STATE.STK.r(:,3)=stkFindData(stkData{1}, 'z ' ) ; 

STATE.STK.V(:,1)=stkFindData(stkData{l}, 'vx' ); 

STATE.STK.V(:,2)=stkFindData(stkData{1}, 'vy' ); 

STATE.STK.V(:,3)=stkFindData(stkData{l}, 'vz' ); 

STATE.STK.ECI=[STATE.STK.r, STATE.STK.v]; 


with the variables, prefixed with STATE.STK, arbitrary assigned. The STK ephemeris 
data is in column format with each row representing a numerical integration step. 
Although not the focus of this research, the user can also pass spacecraft attitude data 
back to the MATLAB Workspace. Eor completeness, format of the attitude code is as 
follows: 
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[stkData,stkNames]=stkReport( '*/Satellite/Target' , 'Attitude 
Quaternions ' ); 

STATE.STK.Qbn(:,iqi + 1)=stkFindData(stkData{1}, 'ql' ); 

STATE.STK.Qbn(:,iqi+2)=stkFindData(stkData{l}, 'q2' ); 

STATE.STK.Qbn(:,iqi + 3)=stkFindData(stkData{l}, 'q3' ); 

STATE.STK.Qbn(:,iqi+4)=stkFindData(stkData{1}, 'q4' ); 

In this STK quaternions format, the first three elements are the vector components and 
the fourth element of the quaternion is the scalar (rotational) term. 

2. Spacecraft Model Validation Results 

Once the independently generated STK spacecraft data is passed into MATLAB, 
the propagated data can be compared against that of the MATLAB/Simulink model. 
Sample average differences between equivalent MATLAB and STK propagation 
methods are reported in Table 10.2. The results of this particular comparison appear to 
validate the MATLAB/Simulink numerical model. These averages were estimated over 
several LEO orbital propagations of 10 spacecraft at various initial conditions. Each orbit 
was propagated for 30 minutes and one minute durations with fixed one second step size. 
The step size introduces floating point errors in the calculation, such that relatively large 
integration step sizes will result in less precision. Also, as the propagation time duration 
increases the error between the propagations slowly and steadily increases. Eor instance, 
if the iteration step size is shortened to 0.5 seconds, the total perturbation position and 
velocity error for a 30 minute duration improves to 1.251 m and 1.895x10'^ m/s, 
respectively. The 30 minute time period was chosen in our simulations as the sample 
duration for a close proximity maneuver. The one minute duration represents an 
estimated time for slow communication or sensor update rates. Errors are a result of 
slight differences in the numerical integration techniques used for orbit propagations and 
numeric precision of constants/coefficients used during the calculation of perturbation 
forces. Overall, the relatively small average error in position and velocity validates the 
MATLAB model. Either the STK HPOP settings or the MATEAB/Simulink model can 
be refined and adjusted in order to achieve equivalent results. Our goal was not to 
achieve identical results, but to validate our spacecraft model by showing similar 
performance. The performance comparison thresholds may vary depending on the user’s 

specific application. Also, some of the computational discrepancies between generated 
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data have periodic characteristics which result in variations of average errors due to the 
time period, or duration, selected. For instance, the average differences between two sets 
of ephemeris over one half of an orbital period may be less then those generated over one 
quarter of an orbital period. 


Perturbation Model 

Average 

Difference 

30 minutes 

1 minute 

Simple Two Body 

Position 

8.237x10-“’ m 

3.391x10"’'’ m 

Velocity 

6.890xlO'*' m/s 

6.082x10"’' m/s 

J2 Perturbation 

Position 

4.650x10"’ m 

1.102x10"' m 

Velocity 

6.908x10"" m/s 

5.346x10"' m/s 

J4 Perturbation 

Position 

8.286x10"’ m 

1.452x10"' m 

Velocity 

1.256x10-' m/s 

7.183x10"' m/s 

Aero Drag 

Position 

1.725x10"'’ m 

1.301x10"' m 

Velocity 

2.904 xlO"’’ m/s 

6.291x10"’'’ m/s 

Solar Drag 

Position 

2.075x10"' m 

6.656x10"' m 

Velocity 

4.267x10"'’ m/s 

3.277xlO"*’ m/s 

Third Body Effects 

Position 

5.721x10"' m 

6.330x10"' m 

Velocity 

1.065x10"" m/s 

3.777x10"® m/s 

Total Perturbation 

Position 

1.360 m 

3.739x10"' m 

Velocity 

2.047x10"' m/s 

3.357x10"" m/s 


Table 10.2 Average Difference of MATLAB and STK Propagation. 


Due to the numerical precision and numerous environmental model subtle 
differences, comparison between different propagations techniques can be challenging. 
The numerical integration and step size precision, mentioned above, is a common source 
of discrepancy. However, orbital propagation algorithms can also vary greatly due to 
small precision and significant digit variations on constants used in near space orbital 
environment. For instance, the most common value for the Gravitational Parameter of 
Earth is // = 398600.4418xlO’m^ [32], however some common applications are based 
on models which use a default value of ju = 398600.4415 xlO^m^ [36]. This appears 
at first glance to be a small change in the last digit, however due to the magnitude of the 
value it can cause large variations in the gravitational effects on the spacecraft model. 
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The orbital discrepancy can be as large as several kilometers, over an orbital duration. 
These challenges are only compounded as the numerous perturbation effects are included 
into the spacecraft dynamics model. Large relative scale differences in computational 
constants can result in loss of desired numerical precision. The spacecraft models used in 
this research is scaled for precision on the order of meters. 

Nevertheless, the relatively short duration for the considered maneuver limits the 
divergence due to perturbations. Minor propagation differences can cause difficulty for 
optimization based algorithms that depend on long duration path propagation and 
tracking. However, close proximity feedback based spacecraft control does not usually 
see these large perturbation changes. The difference between propagators for the 
duration of a sensor measurement cycle, typically less then one minute, are much smaller 
then those above. The short duration of the close proximity maneuvers and the 
disturbance rejection of feedback control make these systems rather tolerant of orbital 
perturbations [63] [64]. Onboard sensors are capable of providing relative position and 
velocity data at rates much less then one minute. Therefore, the values in right column of 
Table 10.2 are well within an acceptable range for the close proximity maneuvers. 
Position error for short durations is approximately two millimeters. This propagation 
position error may be a good metric for determining the size requirements of on-orbit 
docking mechanisms. However, if onboard sensors and communication fail the 
propagation model used by each spacecraft should be standardized so that position 
estimation can be maintained during such blackout periods 

In the future, the multiple spacecraft model can be compared against other high 
fidelity six DOF spacecraft models. For nano-satellite comparison, the Autonomous 
Rendezvous and Rapid Turnout Experiment Maneuverable Inspection Satellite 
(ARTIMIS) testbed at the University of Texas at Austin may be a likely candidate [78]. 
However, determining the cause of variations in model performance is very labor 
intensive. Simulation comparison is especially difficult if different numerical integration 
methods, synchronization, and software or coding are implemented. Ultimately, some 
variation in simulation performance is expected as perturbations are modeled based on 
empirical estimates of data. Precise model refinement would require feedback of 
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operational or test data for actual spacecraft, or components, in the same size and mission 
range. This consideration is not necessary or practical at the algorithm developmental 
level. 


3. Spacecraft Model Propagation Challenges 

Due to the numerical precision and numerous environmental model subtle 
differences, comparison between different propagations techniques can be challenging. 
The numerical integration and step size precision, mentioned above, is a common source 
of discrepancy. However, orbital propagation algorithms can also vary greatly due to 
small precision and significant digit variations on constants used in near space orbital 
environment. For instance, the most common value for the Gravitational Parameter of 
Earth is // = 398600.4418x 10‘^m^ [32], however some common applications are 

based on models which use a default value of //= 398600.4415 xlO'^m^ [36]. This 
appears at first glance to be a small change in the last digit, however due to the magnitude 
of the value it can cause large variations in the gravitational effects on the spacecraft 
model. The orbital discrepancy can be as large as several kilometers, over an orbital 
duration. These challenges are only compounded as the numerous perturbation effects 
are included into the spacecraft dynamics model. Large relative scale differences in 
computational constants can result in loss of desired numerical precision. The spacecraft 
models used in this research is scaled for precision on the order of meters, not the 
standard orbital scale of kilometers. 

The advantage of close proximity research is in the relatively short duration of the 
spacecraft maneuvers. This relatively short duration, usually on the order of a quarter of 
an orbit, limits the divergence due to perturbations. Minor propagation differences can 
cause difficulty for optimization based algorithms that depend on long duration path 
propagation and tracking. However, feedback based control does not see these large 
perturbation changes. The difference between propagators for the duration of a sensor 
measurement cycle, typically less than one minute, are much smaller then those above. 
The short duration of the close proximity maneuvers and the feedback of relative position 
information make the APE controller rather tolerant to perturbation forces. Onboard 
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sensors are expected to provide relative position and velocity data at rates much less than 
one minute. Therefore, the values in right column of Table 10.2 are well within an 
acceptable range for the close proximity maneuvers. Position error for short durations is 
approximately two millimeters. This propagation position error may be a good metric for 
determining the size requirements of on-orbit docking mechanisms. However, if onboard 
sensors and communication fail the propagation model used by each spacecraft should be 
standardized so that position estimation can be maintained during such blackout periods. 


C. SPACECRAFT MODEL ANIMATION 

Accurate 3D simulation and visualization of multiple spacecraft close proximity 
maneuvers support engineering analysis; communicating the complexity and risk in these 
operations [76]. Multiple spacecraft maneuvers were initiated via MATLAB and 
animated via STK. The MATLAB initiates an STK TCP/IP connection to send 
commands over a specified port. The STK scenario is developed by executing STK 
commands. Spacecraft dimensional models for STK visualization must be pre-written in 
STK Modeler. Once the Model is loaded and the scenario is established data can be 
passed between STK and MATLAB. The coding for versatile spacecraft model selection 
and visualization can be extensive and complex. For this research visualization was 
limited to basic cubic, spherical, and cylindrical spacecraft models. More detailed 
models are available for use and reference in the STK model library. 

Typically all control algorithm development is performed in MATLAB with STK 
used for performance visualize and verification. Comparison of control algorithm 
metrics is can be numerically evaluated in MATLAB. However, even with extensive 
metrics and spacecraft state information plots, standardized visualization of the 3D 
spacecraft environment is desired. In order to animate the maneuver, the MATLAB 
generated spacecraft ephemeris can be formatted and passed to STK to be viewed in 2D 
or 3D. The final STK scenario is very useful for visualizing the multiple spacecraft 
maneuver and assist in understanding and troubleshooting any issues that may arise. A 
snapshot of a sample 3D animation of multiple spacecraft is shown in Figure 10.4. Three 
cubic spacecraft are converging toward a common Target spacecraft, with a spherical 

obstacle in the background on the right. Both the view point and proximity are easily 
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adjusted with a standard computer mouse click and drag interface. The animation time 
steps can also be easily adjusted via toolbar buttons. These features allow the user to 
continually evaluate multiple spacecraft simultaneously performing proximity 
maneuvers. For instance as spacecraft are maneuvering through free-space the time step 
can be kept relatively high, but can be slowed down as spacecraft converge. Also, the 
viewpoint may be changed as one spacecraft passes close to an obstacle or eclipses 
another spacecraft. This is extremely useful for evaluating simultaneous docking 
maneuvers of multiple spacecraft. Finally, the user can zoom in on the any desired 
spacecraft to visually check that maneuver constraints, such as spatial safety margins, are 
properly maintained. 



Figure 10.4 Sample 3D STK View of Animation Frame. 

1. STK Spacecraft Model 

For standard STK model assignment, a STK .mdl file must be available. This file 
is assigned to the desired STK satellite object using the following command; 
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stkExec(conID, ['VO */Satellite/Target Model Filename 

"C:\Userflies\STK\Scenario\mdl_cube' ,num2str(isi), '.mdl"'] ); 

where isi is a numerical variable which was used to distinguish multiple uses of the 
same basic model file. This is useful for the modeling of multiple spacecraft with the 
same basic model properties. A sample spacecraft model filename may be 
mdl_cubel .mdl. This would assign the STK satellite object named, Target, with the 
model spacecraft model described by mdl_cubel. mdl. The VO command can be used 
to call a general spacecraft model, which is available in STK. The spacecraft model can 
be selected from a default list, located at C:\Program Files\AGI\STK 
7\STKData\VO\Models\Space, or a custom user generated spacecraft graphical 
model. These models can be manipulated an explored by opening the STK Modeler 
application and loading the desired spacecraft model. Some spacecraft advanced 
visualization setting can be modified from the basic model, by executing selected VO 
commands. These modifications may include scaling, labeling, and sensor descriptions. 


a. STK Graphical Model Development 

Three simple and distinctive spacecraft models were developed in this 
research. They are based on a simple spherical, cylindrical, or cubic spacecraft designs. 
Once again the fprintf command serves as our primary formatting tool. In this 
example the cubic shape serves as the primary body component with supplemental 
docking ports and thrusters components. The sample cubic spacecraft model code is 
presented for discussion and reference as follows: 


filename=[ 'C:\Userflies\STK\Scenario\mdl_cube' ,num2str(isi), '.mdl' ] ; 

mdl_cube=fopen(filename, 'wt'); 

fprintf(mdl_cube, 'Component faceXn'); 

fprintf(mdl_cube, '\t Polygon\n' ); 

fprintf(mdl_cube, '\t\t SpecularityXt 0.2\n'); 

fprintf(mdl_cube, '\t\t ShininessXt 25.6Xn'); 

fprintf(mdl_cube, 'XtXt TranslucencyXt 0.025Xn'); 

fprintf(mdl_cube, 'XtXt FaceColor burlywoodXn' ); 

fprintf(mdl_cube, 'XtXt BackfaceCullableXt NoXn'); 

fprintf(mdl_cube, 'XtXt Translucency 0.025Xn'); 

fprintf(mdl_cube, 'XtXt NumVertsXt 4Xn'); 

fprintf(mdl_cube, 'XtXt DataXn'); 

fprintf(mdl_cube, ['XtXt -' ,num2str(L/2), 'Xt ' ,num2str(L/2), 'Xt - 
' ,num2str (L/2),'Xn']); 
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fprintf(mdl_cube,[ '\t\t ' ,num2str(L/2), '\t ' ,num2str(L/2), '\t - 
' ,num2str (L/2),'\n']); 

fprintf(mdl_cube, ['\t\t ' ,num2str(L/2), '\t ' ,num2str(L/2), '\t 
' ,num2str (L/2),'\n']); 

fprintf(mdl_cube, ['\t\t -' ,num2str(L/2), '\t ' ,num2str(L/2), '\t 

' ,num2str (L/2),'\n']); 

fprintf(mdl_cube, '\t EndPolygon\n' ); 

fprintf(mdl_cube, 'EndComponentVn' ); 

fprintf(mdl_cube, 'Component leftXn'); 
fprintf(mdl_cube, '\t Rotate 0 0 90\n'); 
fprintf(mdl_cube, '\t Refer\n' ); 
fprintf(mdl_cube, '\t Component face\n'); 
fprintf(mdl_cube, '\t EndRefer\n' ); 
fprintf(mdl_cube, 'EndComponentXn' ); 

fprintf(mdl_cube, 'Component topin'); 
fprintf(mdl_cube, '\t Rotate -90 0 0\n'); 

fprintf(mdl_cube, ['\t Translate 0 0 ' ,num2str(L), '\n'] ); 

fprintf(mdl_cube, '\t Refer\n' ); 

fprintf(mdl_cube, '\t Component faceXn'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'EndComponentXn' ); 

fprintf(mdl_cube, 'Component backin') ; 

fprintf(mdl_cube,[ 'Xt Translate 0 -' ,num2str(L), ' OXn' ] ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'Xt Component faceXn'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'EndComponentXn' ); 

fprintf(mdl_cube, 'Component rightin') ; 
fprintf(mdl_cube, 'Xt Rotate 0 0 90Xn'); 

fprintf(mdl_cube, [ 'Xt Translate ' ,num2str(L), ' 0 OXn ' ] ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'Xt Component faceXn'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'EndComponentXn' ); 

fprintf(mdl_cube, 'Component bottomin') ; 
fprintf(mdl_cube, 'Xt Rotate -90 0 OXn'); 
fprintf(mdl_cube, 'Xt ReferXn' ); 
fprintf(mdl_cube, 'Xt Component faceXn'); 
fprintf(mdl_cube, 'Xt EndReferXn' ); 
fprintf(mdl_cube, 'EndComponentXn' ); 


This section of code generates a primary box shape by defining and rotating a simple 
face. The variable, L, represents user defined dimensional parameters which can be 
implemented in the model. In this code the face component was defined first. Next, this 
component was translated and rotated as necessary to determine all size sides of the 
cubic. 
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Additional components can be defined, such as a simple cone. The code 
for a conic docking component is as follows: 


fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube. 


'Component dockingpointVn' ); 

'\t Cylinder\n ') ; 

'\t\t TranslucencyXt 0.25\n'); 

'\t\t BackfaceCullableXt Yes\n'); 

'\t\t NumSidesVt 32\n'); 

['\t\t FacelPositionVt ' ,num2str(L/10), ' 0 0\n']); 
['\t\t FacelRadius\t ' ,num2str(L/20), '\n' ]); 

'\t\t FacelNormalVt 1 0 0\n'); 

['\t\t Face2Position\t -' ,num2str(L/10), ' 0 0\n']); 
['\t\t Face2Radius\t ' ,num2str(L/lOO) , '\n'] ) ; 

'\t\t Face2Normal\t -1 0 0\n'); 

'\t EndCylinderVn' ); 

'EndComponent\n\n' ); 


More detailed components can be generated by combining components. 
For instance, a conic thruster can be combined with a flame, as follows: 


fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube. 


'Component Thrust_Cone\n' ); 

'\t Cylinder\n ') ; 

'\t\t Translucency 0.05\n'); 

'\t\t BackfaceCullableXt NoXn'); 

'XtXtFaceColorXt graylGXn'); 

'XtXt NumSidesXt 32Xn'); 

['XtXt FacelPositionXt -' ,num2str(L/10), ' 0 OXn' 
['XtXt FacelRadiusXt ' ,num2str(L/40), 'Xn' ]); 
'XtXt FacelNormalXt -1 0 OXn'); 

['XtXt Face2PositionXt ' ,num2str(L/10), ' 0 OXn'] 
['XtXt Face2RadiusXt ' ,num2str(L/200), 'Xn'] ); 
'XtXt Face2NormalXt 1 0 OXn'); 

'Xt EndCylinderXn' ); 

'EndComponent Xn' ); 


fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube. 


'Component Thrust_FlameXn' ); 

'PolygonMeshXn' ); 

'Xt FaceColorXt blackXn'); 

'Xt FaceEmissionColorXt graylOOXn'); 

'Xt NoDiffuseLightingXn' ); 

'Xt SmoothShadingXt NoXn'); 

'Xt TranslucencyXt OXn'); 

'Xt SpecularityXt OXn'); 

'Xt Shininess 5lXn'); 

'Xt TextureXn' ); 

'XtXt RGB flametex-whiteXn' ); 

'XtXt Alpha flamealpha3Xn' ); 

'XtXt Farm AAXn'); 

'Xt EndTextureXn' ); 

'Xt NumVerts 12Xn'); 

'Xt DataTxXn' ); 

'XtXt OXt O.OOlXt -l.SlSXt OXt OXn'); 
'XtXt OXt O.OOlXt 1.287Xt OXt iXn'); 
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fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube, 
fprintf(mdl_cube. 


'\t\t -10.55\t -O.OOlXt 1.287\t l\t l\n'); 
'\t\t -10.55\t -O.OOlXt -1.518Xt iXt OXn'); 
'XtXt OXt 1.2lXt 0.586Xt OXt OXn'); 

'XtXt OXt -1.21Xt -0.817Xt OXt iXn'); 

'XtXt -lO.SSXt -1.2lXt -0.817Xt iXt iXn'); 
'XtXt -lO.SSXt 1.2lXt 0.586Xt iXt OXn'); 
'XtXt OXt -1.2lXt 0.586Xt OXt OXn'); 

'XtXt OXt 1.2lXt -0.817Xt OXt IXn'); 

'XtXt -lO.SSXt 1.2lXt -0.817Xt iXt iXn'); 
'XtXt -lO.SSXt -1.2lXt O.S86Xt iXt OXn'); 

'Xt NumPolys 3Xn'); 

'Xt PolysXn' ); 

'4321 OXn' ); 

'4 7 6 S 4Xn' ); 

'4 11 10 9 8Xn' ); 

'Xt EndPolygonMeshXn' ); 

'EndComponent Xn' ); 


The Thrust_Flame component uses textures, such as flametex-white and 
f lamealphaS , to render a more complex image. For detailed discussion of textures refer 
to AGI [36]. 


Once the basic components are defined, they can be reproduced within a loop. 
In this example six thruster cones, based on Thrust_Cone, will be placed on the center of 
each face of the cubic spacecraft. Their respective rotations and positions can be determined 
from a user defined matrix, such as thrustr (thrusti ,: ) and 
thrustl(thrusti, :) . 


for thrusti=l:6 %thruster body positions (-z,-x,x,z 

fprintf(mdl_cube,[ 'Xt Component 
Thrust_Cone' ,num2str(thrusti), 'Xn' ]); 

fprintf(mdl_cube, ['XtXt RotateXt ', thrustr(thrusti,:), 
fprintf(mdl_cube, ['XtXt TranslateXt ', thrusti(thrusti, 
fprintf(mdl_cube, 'XtXt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component Thrust_ConeXn' ); 
fprintf(mdl_cube, 'XtXt EndReferXn' ); 
fprintf(mdl_cube, 'Xt EndComponentXn' ); 


-y.y) 


Xn' ]); 

), 'Xn']); 


fprintf 

num2str 

fprintf 

fprintf 

fprintf 

fprintf 

fprintf 

fprintf 

num2str 

fprintf 

fprintf 

fprintf 

end 


(mdl_cube, 
(thrusti), 
(mdl_cube, 
(mdl_cube, 
(mdl_cube, 
(mdl_cube, 
(mdl_cube, 
(mdl_cube, 
(thrusti), 
(mdl_cube, 
(mdl_cube, 
(mdl_cube. 


[ 'Xt Component Thrust_Elame' , 

'Xn'] ); 

['XtXt RotateXt ', thrustr(thrusti, :),' Xn ']) ; 
['XtXt TranslateXt ', thrusti(thrusti, :),' Xn ']) ; 
'XtXt ReferXn' ); 

'XtXt Component Thrust_FlameXn' ); 

'XtXt EndReferXn'); 

['XtXt ArticulationXt Thrust_Flame' , 

'Xn'] ); 

'XtXt uniformScaleXt SizeXt OXt OXt iXn'); 

'XtXt EndArticulationXn' ); 

'Xt EndComponentXn'); 
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The STK Articulation command is further discussed in Chapter X.C.2.b. 

Now that the primary spacecraft components, including the cubic faces, 
docking ports, and thrusters, are defined, they can be assembled. Each of the six cubic 
components will make the basic spacecraft shape. Next, a thruster with flames 
articulation will be centered on each face of the cubic spacecraft. Finally, the desired 
number of docking ports will be added based on the total number of spacecraft. For 
instance, the Target spacecraft, numerically represented by isi=l will have six docking 
ports. While, all other Chase spacecraft will have single color coordinated docking port 
assigned from user defined matrices, such as dockc (isi-1, :), 
dockr (dockpi, :), and dockl (dockpi , :), for color, rotation, and position, 
respectively. Each component is referred onto the graphical spacecraft model as follows: 


fprintf(mdl_cube, 'Component spacecraftVn' ); 

fprintf(mdl_cube, 'Root\n' ); 

fprintf(mdl_cube, '\t Refer\n' ); 

fprintf(mdl_cube, '\t\t Component faceXn'); 

fprintf(mdl_cube, '\t EndRefer\n' ); 

fprintf(mdl_cube, '\t Refer\n' ); 

fprintf(mdl_cube, '\t\t Component leftXn'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component topXn'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component backin'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component rightin'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component bottomin'); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

for thrusti=l:6 %% thruster positions (-z,-x,x,z,-y,y) 
fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, ['XtXt Component Thrust_Cone' ,num2str(thrusti), 
'in'] ); 

fprintf(mdl_cube, 'Xt EndReferXn' ); 
fprintf(mdl_cube, ' Xt ReferXn' ); 

fprintf(mdl_cube, ['XtXt Component Thrust_Flame ', num2str(thrusti), 
'in'] ); 

fprintf(mdl_cube, ' Xt EndReferXn' ); 
end 

if isi==l 

for dockpi=l:6 
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fprintf(mdl_cube, '\t Refer\n' ); 

fprintf(mdl_cube, '\t\t Component dockingpointXn' ); 

fprintf(mdl_cube,[ '\t\t FaceColorXt ' ,dockc(dockpi,:), 'Xn' ]) 

fprintf(mdl_cube,[ 'XtXt RotateXt ' ,dockr(dockpi,:), 'Xn'] ); 

fprintf(mdl_cube,[ 'XtXt TranslateXt ' ,dockl(dockpi,:), 'Xn' ]) 

fprintf(mdl_cube, 'Xt EndReferXn' ); 

end 

elseif isi<=7 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component dockingpointXn'); 
fprintf(mdl_cube, [ 'XtXt FaceColorXt ' ,dockc(isi-1, :), 'Xn ' ]); 
fprintf(mdl_cube, [ ' XtXt RotateXt ' ,dockr(3, :), 'Xn' ]); 
fprintf(mdl_cube, [ 'XtXt TranslateXt ' ,dockl(3, :), 'Xn ' ]); 
fprintf(mdl_cube, 'Xt EndReferXn' ); 
else 

fprintf(mdl_cube, 'Xt ReferXn' ); 

fprintf(mdl_cube, 'XtXt Component dockingpointXn'); 
fprintf(mdl_cube,[ 'XtXt FaceColorXt ' ,dockc(1,:), 'Xn' ]); 
fprintf(mdl_cube,[ 'XtXt RotateXt ' ,dockr(3,:), 'Xn' ]); 
fprintf(mdl_cube,[ 'XtXt TranslateXt ' ,dockl(3,:), 'Xn'] ); 
fprintf(mdl_cube, 'Xt EndReferXn' ); 
end 

fprintf(mdl_cube, 'EndComponentXnXn' ); 
fclose(mdl_cube); 


The entire graphical modeling section of code can be called as one 
MATLAB file, such as model_cubic_dock_thrust .m. This allows this graphical 
modeling to be isolated from the other MATLAB-STK simulation interface code. 


b. STK 3D Visualization Options 

Several graphics and VO parameters can be tailored for clear visualization 
of multiple spacecraft maneuvers. The primary scaling of the model view is 

accomplished by selecting the following VO command; 

stkExec(conID, 'VO */Satellite/Target ScaleModel Ratio'); 

where Ratio is the number of times larger that the spacecraft model appears during 
animation. For instance, spacecraft being observed from a lunar perspective may nee to 
appear hundreds of times larger in order to appear form this perspective. In close 
proximity operations, Ratio=l was selected in order that spacecraft appear strictly as 
they are modeled. The label and graphic features of the spacecraft can be modified as 
follows: 
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stkExec (conID, 'VO */Satellite/Target LabelOffsetInPixels Off'); 
stkExec(conID,[ 'VO */Satellite/Target LabelXYZ 0.0 ' ,num2str(L), ' 
0 . 0 '] ) ; 

StkExec(conID, 'Graphics */Satellite/Target Basic Color cyan LineWidth 

4.2') ; 


The orbital path is typically shown for spacecraft propagation. However, for multiple 
spacecraft these projections can result in visual confusion. These orbital pass projections 
can be modified as follows: 


StkExec(conID, 'VO */Satellite/Target Pass OrbitLead None'); 

StkExec(conID, 'VO */Satellite/Target Pass OrbitTrail None'); 

In addition to labeling the spacecraft, inclusion of a body axis may be useful for attitude 
reference. A 3D spacecraft body reference frame can be added to a model as follows: 

StkExec(conID, 'VO */Satellite/Target SetVectorGeometry Data 
FixVectorAxesScaleToModel On Scale 0.01'); 

StkExec(conID, 'VO */Satellite/Target SetVectorGeometry Modify 
"Satellite/Target Body Axes" Show On ShowLabel On Color Gold 
ArrowType 3D'); 


Besides straightforward labeling and graphics modifications, additional 
sensors may be desired. These sensors may represent actual ranging and communication 
devices, or useful visual projections for engineering analysis. Each sensor must be 
uniquely named, so that STK dos not confuse identically named sensors in the same 
scenario. In this research a conical ranging and docking sensor was added to each Chase 
spacecraft. The ranging sensor represents the local sensor view from the Chase 
spacecraft to the Target spacecraft, or position. The simple conic range sensor is 
assigned as follows: 


StkExec(conID, ['Location */Satellite/', Chase, '/Sensor/Range' ,Chase, ' 
Fixed Cartesian ' ,num2str(L/2+0.001), ' 0 0']); 

StkExec(conID,[ 'Define */Satellite/' ,Chase, '/Sensor/Range' ,Chase, ' 
SimpleCone 0.25']); 

StkExec(conID, ['Point ^/Satellite/', Chase, '/Sensor/Range' ,Chase, ' 
Targeted Tracking Satellite/Target Hold']); 

StkExec(conID,[ 'VO */Satellite/' ,Chase, '/Sensor/Range' ,Chase, ' 
Projection SpaceProjection ' ,num2str(rsw_norm(1,isi+1)-L)]); 

StkExec(conID,[ 'Graphics */Satellite/' ,Chase, '/Sensor/Range' ,Chase, ' 
SetColor cyan']); 


198 


First, the sensor is positioned and defined on a spacecraft with an exclusive name. In this 
example, Chase is an alpha-numeric label unique to each Chase spacecraft. Next, the 
pointing direction of the sensor is determined. The range sensor continually points 
toward the Target spacecraft and adjusts its projection based on the relative distance from 
the Target, represented by rsw_norm. This modification of the projection may not be 
necessary for most sensor applications. Finally, the color of the sensor projection can be 
modified as desired. Similarly, a docking sensor may be added to the spacecraft 
graphical model as follows; 


stkExec(conID, ['Location */Satellite/', Chase, '/Sensor/Dock' ,Chase, ' 
Fixed Cartesian ' ,numCstr(L/2+0.001), ' 0 0']); 

StkExec(conID, ['Define */Satellite/ ', Chase, '/Sensor/Dock ', Chase, ' 
SimpleCone 45.0']); 

StkExec(conID, ['Point */Satellite/', Chase, '/Sensor/Dock' ,Chase, ' 
Fixed YPR 321 0 90 0' ]); 

StkExec(conID,[ 'VO */Satellite/' ,Chase, '/Sensor/Dock' ,Chase, ' 
Projection SpaceProjection ' ,num2str(L/2)]); 

StkExec(conID,[ 'Graphics */Satellite/' ,Chase, '/Sensor/Dock' ,Chase, ' 
SetColor ' ,dockc(isi-1,:)]); 


The docking sensor is located at the same position, but has a much wider cone of 45 
degrees. The docking sensor points in a fixed direction from the spacecraft body axis, 
with a limited projection of half of the spacecraft’s length. Next, the sensor projection 
color is selected from a matrix, called dockc. These simple sensors are useful in 
visually representing ranges and FOV as spacecraft are animated through their dynamics 
and kinematics. A sample cubic spacecraft with labels and sensor projections is shown in 
Figure 10.5. The three body axis are shown in yellow with labels. The green projection 
on the right side of the spacecraft is the docking sensor projection and the small cyan line 
extending to the bottom left is the range sensor. The white protrusions on ht top and right 
are thruster firings. 
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One additional sensor projeetion was useful in visualizing spaeeeraft 
exelusion regions during elose proximity operations. A boundary sensor projection was 
placed around regions that spacecraft were intended to avoid. These areas were referred 
to as obstacles and labeled with Obstcl. Spherical projections encompassing these 
regions can be generated using conical sensor defined with 360 degrees FOV, such as 
follows: 


stkExec(conID, ['Define */Sat ellite/', Obstcl, '/Sensor/Bound ', Obstcl, ' 
Conical 0 180 0 360']); 

StkExec(conID, [ ' SetConstraint 

*/Satellite/ ', Obstcl, '/Sensor/Bound' ,Obstcl, ' Range Max 
' ,num2str(OBS.DoL)]); 

StkExec(conID, ['Location 

*/Sateilite/ ', Obstcl, '/Sensor/Bound' ,Obstcl, ' Fixed Cartesian 0 0 
O'] ) ; 

StkExec(conID,[ 'VO */Sat ellite/', Obstcl, '/Sensor/Bound' ,Obstcl, ' 
Projection SpaceProjection ' ,numCstr(OBS.DoL)]); 

StkExec(conID, ['Graphics 

^/Satellite/ ', Obstcl, '/Sensor/Bound' ,Obstcl, ' SetColor red']); 
StkExec(conID,[ 'VO */Sat ellite/', Obstcl, '/Sensor/Bound' ,Obstcl, ' 
Translucency 95']); 

StkExec(conID,[ 'VO */Sat ellite/', Obstcl, '/Sensor/Bound' ,Obstcl, ' 
TranslucentLines On']); 
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The spherical sensor projections were made translucent in order not to 
obscure the view of maneuvering spacecraft. A sample spherical boundary sensor 
projection about a spherical object is shown in Figure 10.6. The sensor projection is in 
translucent red with the Chase spacecraft in the upper right. These sensors allowed for 
robust engineering evaluation of collision avoidance function of spacecraft control 
algorithms. From these basic examples, the user can develop a vast array of 
representative sensors. Sensor projections ranges and colors can be varied based on 
logical condition. For instance, ranging sensors may be programmed to modify their 
color as the distance to a target varies. This may allow for quick visual analysis of 
multiple spacecraft maneuver performance. However, users may also want to terminate 
some projections during presentations, since multiple overlapping sensors may be 
visually confusing to some audiences. Ultimately, the animation graphic and level of 
detail can be tailored from the core MATLAB-STK simulation interface code described 
in this research. 



Figure 10.6 Sample Spherical Boundary Sensor. 
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2. Formatting MATLAB Data for STK Files 

The visualization and animation of the MATLAB/Simulink data with the desired 
spaeecraft model is a useful tool for engineering analysis and evaluation. In partieular, it 
was vital in the development of a multiple spaeeeraft eontrol algorithm during 
simultaneous close proximity operations [63][64], The 3D representation for spacecraft 
during collision avoidance maneuvers enabled effective troubleshooting. For instance, 
the undesired clipping of the comers of Chase spacecraft while converging to the Target 
spacecraft during docking maneuvers was easily shown in STK animation. Due to this 
evaluation, the control algorithm collision avoidance logic was modified and its 
performance robustness was improved. Additionally, visualization of the simultaneous 
multiple spacecraft maneuvers in STK allows for logical point of view changes and 
animation rate changes. These capabilities allow for clear presentation of spacecraft 
dynamics and control results to those interested multiple spacecraft operations. In fact, 
STK scenarios can be easily edited and transferred into video format for presentations. 
Modifications of camera key frames (camera position and angle) and animation rates 
allow versatile video composition. The video fde’s format, encoding, bit rate, and size 
can all be selected in STK 7.1, or later versions [36]. The basic STK scenario can be 
modified and animated to meet the individual user’s requirements. Ultimately, the same 
STK scenario used for engineering evaluation can serve as a demonstration of 
engineering capabilities to a wider audience. 


a. Ephemeris and Attitude Data 

The animation of the spacecraft can be based on the data generated by 
STK or the MATLAB/Simulink model. STK generated data is already in the proper 
format. However, MATLAB/Simulink data must be proper formatted in order to be used 
by STK. STK allows for several specific data file formats, such as STK ephemeris and 
attitude files. These files can be properly written by executing short MATLAB scripts 
utilizing fprintf commands. A sample STK ephemeris file can be created in 
MATLAB as follows: 
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filename= 'C:\Userflies\STK\Scenario\ephemeris_t.e' ; 

ephemeris=fopen(filename, 'w' ); 

fprintf (ephemeris, 'stk..v.4.3\n\n') ; 

fprintf(ephemeris, 'BEGIN Ephemeris\n\n' ); 

fprintf(ephemeris, 'NumberOfEphemerisPointsXt 

%6.Of\n' ,length(STATE.time)); 

fprintf(ephemeris, 'ScenarioEpochVt %s\n' ,num2str(para_now)); 

fprintf(ephemeris, 'InterpolationMethodXt LagrangeXn' ); 

fprintf(ephemeris, 'InterpolationOrderXt 5Xn' ); 

fprintf(ephemeris, 'CentralBodyXt EarthXn' ); 

fprintf(ephemeris, 'CoordinateSystemXt J2000Xn' ); 

fprintf(ephemeris, 'EphemerisTimePosVelXnXn' ); 

data=[STATE.time(:,l),r(:,ici+l:ici+3),v(:,ici+l:ici+3)] 

fprintf(ephemeris, '%16.14eXt %16.14eXt %16.14eXt %16.14eXt %16.14eXt 

%16.14eXt %16.14eXn' ,data); 

fprintf(ephemeris, 'Xn END EphemerisXn' ); 

fclose(ephemeris); 


The STK ephemeris file name and path, in the first line, are arbitrary, as long as it has the 
proper . e suffix. Although, it is useful to save all related STK seenario files in the same 
folder. The first few fprintf lines are used to establish the desired reference frame and 
interpolation for the data. The variables STATE.time and para_now are used to 
synchronize the data points with the time constraints, as discussed in Chapter I. C. 2. The 
data is formatted into seven columns, represented by time, position vector, and velocity 
vector. Each row represents an iteration, or stepsize, of the data. Similarly, a STK 
quaternion attitude fide can be created as follows: 


filename2= 'C:XUserfliesXSTKXScenarioXattitude_t.a' ; 

attitude=fopen(filename2, 'w' ); 

fprintf(attitude, 'stk.v.S.OXnXn'); 

fprintf(attitude, 'BEGIN AttitudeXnXn' ); 

fprintf(attitude, 'NumberOfAttitudePointsXt 

%6.OfXn' ,length(STATE.time)); 

fprintf(attitude, 'ScenarioEpochXt %sXn' ,num2str(para_now)); 

fprintf(attitude, 'BlockingFactorXt 20Xn' ); 

fprintf(attitude, 'InterpolationOrderXt 5Xn' ); 

fprintf(attitude, 'CentralBodyXt EarthXn' ); 

fprintf(attitude, 'CoordinateAxesXt J2000Xn' ); 

fprintf(attitude, 'AttitudeTimeQuaternionsXnXn' ); 

data=[STATE.time(:,1),Qbn(:,iqi+1:iqi+4)]'; 

fprintf(attitude, '%16.14eXt %16.14eXt %16.14eXt %16.14eXt 

%16.14eXn' ,data); 

fprintf(attitude, 'Xn END AttitudeXn' ); 
fclose(attitude); 
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The STK attitude file name and path are arbitrary, as long as it has the proper . a suffix. 
This sample STK attitude file is formatted into five eolumns with time, the three element 
quaternion veetor, and the quaternion scalar (rotational) term. 

b. STK Model Articulation Files 

The articulation of a STK spacecraft graphic model component is executed 
by padding formatted data to the STK scenario. The articulation file must have the 
. sama suffix, be named the same as the spacecraft object, such as Target, and saved 
in the same folder as the STK graphical model, such as the mdl_cubel. mdl. A 
sample articulation file for each Thrust Flame component is as follows: 

filename3=[ 'C: \Userflies\STK\Scenario\Target.sama' ]; 
articulate=fopen(filenames, 'wt' ); 
fprintf(articulate, 'STARTTIME\t 0\n\n ' ); 
fprintf (articulate, 'DURATION\t %6.Of\n' ,stoptime); 
fprintf(articulate, 'DEADBANDDURATION\t 0.0\n' ); 
fprintf(articulate, 'ACCELDURATION\t 0.0\n' ); 
fprintf(articulate, 'DECELDURATION\t 0.0\n' ); 
fprintf(articulate, 'DUTYCYCLEDELTA\t 0.0\n' ); 
fprintf(articulate, 'PERIOD\t 0.0\n' ); 

fprintf (articulate, 'ARTICULATION\t Thrust_Flame#\n' ); 
fprintf(articulate, 'TRANSFORMATION\t Size\n' ); 
fprintf(articulate, 'STARTVALUE\t 0\n' ); 
fprintf(articulate, 'ENDVALUE\t l\n' ); 
fclose(articulate); 

where Thrust_Flame# is the name of the component being articulated and Size is 
the magnitude of the articulation desired. This strict formatting is limiting due to the start 
and end values of each articulation iteration. A more useful presentation, defines the 
articulation in a spreadsheet which loops through the time increments, ti, and assigns 
the desired articulation magnitude for each thruster, til, on the spacecraft. Sample 
articulation code is as follows: 


for ti=l:length(time) 
for tii=l:6 

filename3=[ 'C: \Userflies\STK\Scenario!Target.sama' ]; 
articulate=fopen(filename3, 'wt' ); 
fprintf (articulate, 'SPREADSHEET\n\n' ); 

fprintf(articulate, [ 'ARTICULATION, ' ,numlstr((ti-1)*stepsize), ' , ' , 
numlstr(stepsize), ',0,0,0,0,0,Thrust_Flame' ,numlstr(tii), ',Size,' , 
numlstr(tmag(ti,tii)), ', ', numlstr(tmag(ti+l,tii)),'\n']); 
end 
end 

fclose(articulate); 
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The spreadsheet data may be based on desired thruster tiring magnitude at eaeh time step 
of the simulation. This data may be generated from the MATLAB/Simulink spacecraft 
model’s control algorithm response. The desired control response can be mapped to the 
appropriate thruster. This mapped response can be assigned into the articulation file to 
animate each spacecraft’s thrusters firing during the multiple spacecraft maneuver. 

It may be useful to include a reload command at the end of the articulation 
file assignment to ensure that all components have the current data. This reload 
command is of the form; 

stkExec(conID, 'VO */Satellite/Target ReloadArticFile' ); 

The articulation of multiple components can be useful in generating complex spacecraft 
models. Articulating components may include such components as thrusters, solar 
panels, and momentum exchange devices. Although visually stimulating, it may be more 
practical to only articulate the necessary components for each engineering application. 
The articulation sections may be commented out when not necessary. Once performance 
is complete, the articulation may be reinstated for presentations. 

D. MATLAB-STK CONCLUSION 

A MATLAB-STK simulation interface was developed for spacecraft model 
validation and visualization. Both MATLAB/Simulink and STK have the capability to 
perform high precision orbital propagation. Utilizing this capability, the spacecraft 
model validation is conducted by comparing various spacecraft state date, such as 
position and velocity. Each modular perturbation component can be independently 
enabled for evaluation. The desired level of model similarity and accuracy can be 
achieved by matching all parameters that are being considered by the designer. Once the 
model is validated, the spacecraft can be animated by passing spacecraft state data to 
STK. The spacecraft characteristics and parameters, defined in MATLAB/Simulink, can 
be formatted and passed into STK. Utilizing STK via MATLAB is a versatile and 
effective method of animating six DOF spacecraft models. The resulting STK animation 
of spacecraft propagations is useful for engineering evaluation and result presentations. 

This interface was utilized during the development of a multiple spacecraft control 
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algorithm for close proximity operations. Model validation gave eonfidenee in the 
performance results and visualization allowed for straightforward evaluation. The 
animation allowed for immediate identifieation of undesirable performanee such that 
modifications and improvements eould be made to the eontrol algorithm. This 
MATLAB-STK simulation interfaee proved to be a useful tool for the development and 
evaluation of spaeeeraft models and eontrol systems. 
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XI. VHIL EVALUATION 


The challenging characteristics of microgravity orbital space make it difficult to 
simulate in the laboratory environment. So, a series of test methods exist for incremental 
evaluation of the control algorithm. First, the control algorithm processing can be 
implemented on distributed independent hardware. This VHIL testbed represents a set of 
simulated spacecraft which conduct the control algorithm based on their local state 
information. Next, a terrestrial based spacecraft robot testbed can be utilized. These 3D 
spacecraft robots can actual physically actuate and perform limited two-axis translational 
and one axis rotational maneuvers. Finally, the control algorithm can be loaded onto 
actual free flying spacecraft in the orbital environment. 

A. VIRTUAL HARDWARE-IN-THE LOOP (VHIL) TESTBED 

Hardware testing for spacecraft related components and systems are restricted due 
to high costs and environmental limitations. In order to perform incremental testing of 
the developed control algorithm, a VHIL testbed was developed at NFS. This test bed 
consists of a distributed set of independent processors linked in a virtual environment. 
The block diagram for the VHIL testbed is shown in Figure 11.1. 
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Each Vehicle sends their incremental state information to the Primary Simulation Controller. 

Each vehicle receives virtual sensor information based on other Vehicle’s state information 
received by the Primary Simulation Controller. 

Figure 11.1 VHIL Block Diagram. 

Each processor represents a stand alone Chase spacecraft which independently 
propagates its maneuver ephemeris based on a validated spacecraft model. Any external 
data concerning the position and velocity of other spacecraft is passed via UDP. This 
data simulates possible sensor derived information. A desktop computer may serve as the 
Target spacecraft. This desktop may also provide processing of the STK visualization for 
the multiple spacecraft maneuvers. The same software used in the original development 
was used in the distributed control algorithm, primarily MATLAB version 7.3.0267 
(R2006b) and STK version 7.0.1 [36][39]. All MATLAB spacecraft dynamics and 
control functions were converted into embedded MATLAB functions for efficient C code 
programming language. The virtual spacecraft were represented by small mobile 
computers: refer to Table 11.1 for details. The multiple spacecraft VHIL testbed is 
shown in Ligure 11.2. This image shows four distributed processors wirelessly 
connected, with a shared monitor, keyboard and mouse. 
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Characteristics 

Name 

LG-P625E 

Processor 

Intel Pentium M 1.7 GHz 

Hard drive 

100 GB 7200 RPM 

Memory 

1 GB DDR400 

DC Power 

12 V 

Size 

Width 

0.155 m 

Depth 

0.255 m 

Height 

0.055 m 

Mass 

<3.0 kg 

Eeatures 

Thermal Control Technology 

Pan less 

Built-in Wireless 

802.1 IB/G 

Slot Loading Disc Drive 

CDRW with DVD 

Interfaces 

USB 2.0, Pirewire 

Operating Systems 

Microsoft Windows 

XP (80 GB) 

Linux 

Gentoo (20 GB) 


Table 11.1 Virtual Spacecraft Computers. 



Figure 11.2 Multiple Spacecraft Testbed. 


The VHIL distribution of the control algorithm ensures that limited state 
information and reasonable data rates are being utilized. Initial assessments of processor, 
sensor, and timing limitations can be conducted. Each processor independently simulates 
the spacecraft environment. Each virtual spacecraft processing rate can be roughly 
synchronized by using Simulink Simulation Pace blocks to standardize the ratio of 
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simulation seconds to internal clock seconds. This serves as a versatile soft-real time 
technique for implementing the VHIL via Simulink [36]. 

The control algorithm executes based on the autonomous spacecraft’s simulated 
sensor information. This data can be passed directly between Chase processors or routed 
via the Target processor. Drop out of one or more Chase spacecraft processors do not 
negatively affect the VHIL simulation of other spacecraft. UDP communication default 
simply holds the last value successfully received [36]. Initial VHIL testing focused on 
control algorithm development and not sensor or communication link limitations. 
Obviously, obstacle and other Chase spacecraft states must be know at some level in 
order permit collision avoidance. The distributed state data for each Chase spacecraft is 
only centralized at the conclusion of simulations for the purpose of evaluation and 
visualization. The VHIL testing paves the way for additional testing on terrestrial and 
orbiting testbeds. 

B. AMPHIS TESTBED 

In order to fully develop the concepts for close proximity operations of multiple 
spacecraft through hardware-in-the-loop testing, it is necessary to have sufficient space 
for maneuvering of the testbed vehicles. A Proximity Operations Simulator Facility was 
been previously developed [79][80] and will be used for any new systems. This facility 
consists of a 4.9 m by 4.3 m wide Epoxy Floor Surface, approximately 21 square meters, 
used as base for the floatation of spacecraft simulators. The use of air pads on the 
simulators reduces the friction to negligible level. Due to an average residual slope angle 
of -2.6-10 Meg for the floating surface, the average value of the residual gravity 

acceleration affecting the dynamics of floating vehicles is - 1.4 10 ^ m*s‘^. This value of 
acceleration, measured by analyzing the free motion of the Chaser spacecraft simulation, 
is two orders of magnitude lower than the nominal amplitude of the acceleration 
fluctuation obtained during the reduced gravity phases of parabolic flights [81]. 

Given the dimensional constraints of the flat epoxy floor, in order to 
simultaneously operate at least three robotic spacecraft a significantly lighter and smaller 
vehicle is required, as compared to the previous AUDUSS spacecraft robot simulator [79] 
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[80]. Additionally, the NPS SRL research and testing of multiple sensor packages 
requires multiple spacecraft robots. The new generation of robotic spacecraft at the SRL 
is referred to as AMPHIS, as shown in Figure 11.3. 



Router 


Dual 

MSGCMG 


iGPS sensor 


LIDAR 


Djal PC104 


3 Li Ion 
Batteriee 


Dual 

Vecto'able 

Thrusters 


3000 psi 
Air Cylinders 


4 Air Pads 


Figure 11.3 AMPFIIS Robotic Spacecraft Simulator at NPS [82]. 


Recent enhancements in the AMPHIS robotic spacecraft simulator enabled a 
significant reduction in size and mass while providing for easy reconfiguration [82][83]; 
refer to Table 11.2. The details on the AMPHIS robot are presented completeness. The 
AMPHIS is continually evolving as research is conducted at NPS. 
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Size 

Eength and Width 

0.30 m 

Height 

0.69 m 

Mass 

37 kg 

Moment of Inertia 

0.75 kg m^ 

Propulsion 

Propellant 

Air 

Equiv Storage Cap @21 MPa (3000 PSI) 

0.002 m^ 

Operating Pressure 

0.827 MPa 

Thrust per Thruster 

0.28 N 

Attitude 

Control 

MSGCMG Max Torque 

0.668 Nm 

MSGCMG Max Angular Momentum 

0.098 Nms 

Electrical & 
Electronic 

Battery Type 

Eithium-Ion 

Storage Cap @ 28 V 

12 Ah 

Computers 

2 PC-104, Pentium III 

Sensors 

Tiber Optic Gyro Bias 

±207hr 

LIDAR 

SICK 360 ° 

IGPS Sensor Accuracy 

<.050 mm 

Accelerometers Bias 

±8.5x10'^ g 

Eloatation 

Propellant 

Air 

Equiv Storage Cap @21 MPA (3000 PSI) 

.002 m" 

Operating Pressure 

0.35 MPa 

Einear Air Bearings Diameter 

32 mm 


Table 11.2 Key Parameters on the AMPHIS Testbed [82][83]. 


It is worth mentioning that other Hardware-in-the-Loop testbeds are being 
developed world wide. All terrestrial spaceeraft testbed has limitations due to the need to 
imitate the orbital environment. The concurrent development of terrestrial testbeds may 
allow for comparison of algorithm performance results. Besides the NPS SRL, another 
promising satellite formation flying testbed is located at the University of Southampton, 
in the United Kingdom [84]. This facility enables control algorithms to be tested on 5 
DOF mechanical mock-up satellite frames, referred to as ground-based satellite frame 
testing (GSFT) [84]. The control algorithm testing relies on simulated sensor, actuator, 
and physical dynamics. The combination of these virtual components and the GSFT 
overlaps for systematic spacecraft environmental simulations. The accuracy of the 
testbed currently aims at centimeter accuracy, but may be improved for higher precision 
[84]. 
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C. SPHERES TESTBED 

On-orbit testing is a tremendous opportunity for development a multiple 
spaeecraft control algorithm. On board the ISS, testing of the multiple spacecraft close 
proximity control algorithm is now possible with the MIT Space Systems Laboratory’s 
SPHERES [25] [26]. The SPHERES are intended to be used as a test bed for formation 
flight and reconfiguration, as well as autonomous rendezvous and docking technologies 
[26]. The first SPHERES satellite reached the ISS in May 2006 and has begun testing 
scenarios. There are currently three SPHERES available for ISS investigation sessions. 

Before on-orbit testing can be conducted, software modification and 
synchronization would have to be conducted. This would include ground testing with 
MIT Space Systems Eaboratory’s team. Eunding and co-operative support for this work 
has been discussed. DARPA, NASA, MIT, NPS and DoD participants have all shown 
interest. Scheduling and timeline concerns have limited the inclusion of this testing in 
this current body of research. 
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XII. CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

The developed LQR/APF multiple spacecraft proximity control algorithm offers 
desirable performance and establishes a baseline for fuel efficiency while maintaining 
collision free operations. Evaluation of the LQR/APF control algorithm proved it to be 
effective in a broad range of multiple spacecraft close proximity maneuvers, including 
rendezvous and docking maneuvers. In the presence of obstacles and other maneuvering 
spacecraft the LQR/APF exhibited a smooth and effective control response which 
avoided actuator saturation. Therefore, the LQR/APF control algorithm shows efficient 
control response with reliable collision avoidance 

Furthermore, Monte-Carlo method analysis of close proximity maneuvers 
confirms that the LQR/APF control algorithm is a practical candidate for multiple 
spacecraft close proximity operations. Estimates of the mean and standard deviation of 
maneuver Q and Av show that average proximity maneuver LQR/APF control Av 
efficiency is better than that of a highly tuned APF control algorithm. The LQR/APF 
showed a distinct control efficiency improvement on a per spacecraft maneuver basis. 
The standard deviation of the LQR/APF control effort is consistently narrower then that 
of the APF controller. This narrow Av standard deviation is valuable to both spacecraft 
designers and mission planners. It allows effective propellant sizing for close proximity 
operations. It also gives operational planners a useful tool for developing and forecasting 
maneuvers. 

In conjunction with the controller research, a MATLAB-STK simulation interface 
was developed for multiple spacecraft model validation and dynamic environment 
visualization. This MATLAB-STK simulation interface was extensively utilized during 
the LQR/APF control algorithm refinement development. Model validation gave 
confidence in performance results and visualization allowed for straightforward 
evaluation. The STK animation allowed for immediate identification of undesirable 
performance such that modifications and improvements could be made to the control 
algorithm. This MATLAB-STK simulation interface is a useful tool for the development 
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and evaluation of spacecraft models and control systems. The STK animation allows for 
clear presentation of developed control algorithm performance to the spacecraft field. 

Finally, the LQR/APF control algorithm was implemented and evaluated in a 
VHIL configuration. The VHIL utilizes independent processors which simulate 
spacecraft and interact as multiple spacecraft. VHIL control algorithm test results were 
consistent with the LQR/APF performance during the evaluation and analysis presented 
in this research. Successful VHIL experimentation paves the way for future hardware 
implementation on both terrestrial based and orbital spacecraft robot testbeds. 

B. FUTURE WORK 

Further investigation could establish the control algorithm robustness by 
conducting additional Monte-Carlo simulations including a higher number of spacecraft, 
heterogeneous spacecraft, noisy sensor, and more restrictive actuator limitations and 
failures. These studies could address the potential scope of application for this control 
algorithm. Elliptical orbital considerations could be included into the dynamics [43] [42] 
[44]. Specific on-orbit applications need to take into account the actual spacecraft sensor 
and actuators in order to determine if the control algorithm’s performance meets mission 
requirements. 

Refinement of the APF control algorithm could be carried out based on a more 
sophisticated Newton method or conjugate gradient search algorithms. The Steepest 
Descent algorithm initially adopted is based on the system measurement and control rates 
being nearly the same. However, the control effort may need to be two or three times 
faster then the measurement cycle. In this case, the quadratic characteristics of the 
potential functions should allow the faster control rate to achieve more efficient iterative 
minimization. In particular the conjugate gradient Fletcher Reeves iteration scheme 
efficiency is due to varying the step magnitude and direction variation while convergence 
is achieved in two iterations. Refer to Appendix A for alternate search algorithm 
discussion. 

Refinement of the LQR control algorithm to include direct incorporation of 
collision avoidance terms into the state and control effort gain matrices could be 
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investigated. If the cost function can be minimized based on collision avoidance, then it 
should be possible to include additional parameters into the LQR gains which are related 
to collision avoidance. The desire would be to add collision avoidance gains without 
adding additional states. The selection of these gains needs to maintain a geometric 
relationship to collision avoidance, while preserving the desired system response. 
Preliminary success was achieved by incorporating obstacle relative velocity components 
into the LQR state gain matrix. In general, collision avoidance was achieved; however 
the tuning of these parameters proved to be very sensitive as the states approached the 
goal. Therefore, docking precision seen by the hybrid algorithm could not be achieved 
while maintaining collision avoidance. More rigorous study of the gain parameters need 
to be conducted before confidence can be given to this incorporated LQR collision 
avoidance method. Refer to Appendix B for preliminary LQR with collision avoidance 
approach and results. 

Specific actuator limitations in both the translational and attitude control could be 
further investigated. Evaluating controller algorithm performance with limited numbers 
or capabilities in actuation may offer interesting practical maneuver limitations. For 
instance two maneuverable thrusters may be able to perform all desired maneuvers. 
However, specific types of thrusters may be constrained due to pluming regions around 
sensitive payloads. Additionally, the coupling of translational and attitude control via 
thruster actuation is an interesting development. 

In this research, initial goal locations and docking ports were commanded, 

however future work may allows for an outer control loop to determine when and where 

each spacecraft is tasked. Higher level control can be implemented based on neural 

networking or dynamic programming multi-vehicle task assignment algorithms [85][86] 

[1]. However, these algorithms may require additional communication between 

spacecraft for efficient tasking decisions, or voting, to be accomplished. This 

communication requirement may lead to more research into wireless sensor networking 

of multiple spacecraft. In such a networking, micro chips with control with modem 

capabilities, such as Maxstream’s small Xbee chips [87], can be used for communication 

and can also support ranging data. The wireless signal has signal strength data, referred 

to as Received Signal Strength Indicator (RSSI) that can be combined with ranging 
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sensors, such as GPS. The signal strength is determined from the decibel magnitude of 
the last packet. With directional antenna, the decibel reading can be directly transferred 
into ranging data. This additional information can improve range precision and 
supplement when other sensor data. An example of wireless mesh networking on 
autonomous vehicles using the Xbee chip is discussed in detail by Bledsoe in [88]. This 
application can be related to the communication of multiples spacecraft robots on a 
terrestrial testbed during proximity operations. 

The multiple spacecraft control algorithm developed in this research showed 
robust disturbance rejection based on deterministic measurements. However, 
comprehensive performance evaluation including process and sensor noise characteristics 
for specific spacecraft configurations may be desired. Any known sensor statistics can be 
used to re-characterize the stochastic process. This uncertainty can be included into the 
control algorithm by modifying the LQR as a LQG regulator, which can incorporate 
process and measurement noise as Gaussian white noises with covariance. Also the use 
of Gaussian based collision avoidance lends itself to direct inclusion of any measurement 
noise in to the obstacle relative range and velocity. The algorithm structure allows for 
convenient inclusion of known or estimated sensor uncertainties. These uncertainties 
may be particularly suitable to the estimation with Kalman Filtering techniques suggested 
by Cristi et al [89]. Large control iteration time steps, of one second, were maintained to 
allow for fusing or filtering of sensor measurements. Even with noisy state 
measurements the probable shortening of control step iterations will allow for acceptable 
algorithm performance. Although, detailed control algorithm performance evaluation 
based on varying levels of filtered, or individual, sensor noise would offer confidence in 
practical and valid safety margins. 

Full evaluation of docking mechanisms and their related performance 

requirements could be undertaken. Various potential docking mechanisms could be 

incorporated in evaluation of the control algorithm [80][81]. For instance, spatial and 

attitude tolerances may allow for less rigorous terminal control. Also passive docking 

mechanism may actually require Chase spacecraft to engage the Target spacecraft with 

some terminal force. These latch mechanism are envisioned for on-orbit assembly, and 

not for autonomous spacecraft with sensitive payloads [81]. However, even simple 
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adhesive docking ports would require some terminal contact force. The trade-off 
between a firm dock and the jarring of the spacecraft could be researched. 

Most importantly, the close proximity control algorithm can be integrated onto 
spacecraft robots with different vehicle and sensor properties. Ground testing could 
include integration onto the NFS AMPHIS in the Spacecraft Servicing and Robotics 
Laboratory. The laboratory hosts the Autonomous Docking and Spacecraft Servicing 
testbed, which may allow for future testing of on-the-ground dynamics with close 
proximity operations. The AMPHIS is built upon previous research on the original 
AUDUSS and the AUDUSS II [79] [80]. Additionally, ground testing may also be 
conducted at MIT’s Space System Laboratory as preparation for on-orbit flight testing of 
the control algorithm on the SPHERES onboard the ISS [25] [26]. The SPHERES 
testbed enables maturation and validation of the control algorithm in the micro gravity 
environment of the ISS as a risk reduction before integration on future spacecraft. This 
will require incorporation of the SPHERES system details into the control algorithm with 
baseline simulations and ground testing. 
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APPENDIX A: OPTIMAL SEARCH ALGORITHMS 


As discussed during the APF development, previous APF control algorithms have 
used the immediate Steepest Descent approaeh to navigation. This is a reasonable 
assumption as to the best eourse for the Chase spacecraft to proceed. However, the 
implementation of iterative searehes in spaeeeraft dynamies yields only a momentary 
current steep gradient. This direetion of approaeh may not be the most efficient path. 
There are alternate ways for selecting both the seareh direetion and length of the step 
size. In partieular, Newton’s method and Conjugate Gradient optimal seareh algorithms 
appear promising for selection of the seareh direetion and step size varies for Spacecraft 
APF based control algorithms However, in order to utilize these seareh methods the 
spacecraft’s sensor update rates may need to be considered in the eontrols systems 
iterative eonvergence algorithm. While the Steepest Descent method has proved 
extremely effective in this researeh, it may not be the best for all multiple vehicle control 
algorithms. Consideration and comparison of other seareh techniques, sueh as Newton’s 
method and Conjugate Gradient, should be made before a final determination can be 
made. 

A. STEEPEST DESCENT 

Whether a potential or cost function is used as the guiding parameter for the 
control algorithm, some discussion of optimal search is warranted. Broadly speaking, 
seareh proeedures for optimization problems may be divided into three main categories: 
calculus-based, enumerative, and heuristie. Calculus-based proeedures use either 
analytieal or numerie models of the solution space as the basis for the seareh. 
Enumerative or “brute force” proeedures seareh systematieally and do not incorporate 
any sophistieated methods to reduce the problem space or refine proposed solutions. 
Heuristie proeedures attempt to improve on the seareh efficieney of enumerative methods 
without direct incorporation of models of the solution spaee, whieh are often unavailable. 

Performanee surfaces are funetions which are designed such that the quantitative 
value is smallest when the desired performanee is obtained. The global minimum point is 
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considered the equilibrium or optimal point [85]. The path whieh the system searehes for 
and approaehes this minimum point is dependent on the optimization algorithm used. 
The iterative seleetion of the seareh direction and step size varies for eaeh type of 
optimization algorithms. While a wide breadth of search algorithms exists, those most 
applieable to this work are the Steepest Deseent, Newton’s method, and eonjugate 
gradient. Given its effieieney in implementation, the Steepest Deseent method was 
employed in this researeh. 

For a performanee function, V (x ), the general iterative seareh algorithm is of the 

form 

V(x,^i)=V(x,+AJ^,) (A.l) 

where x^ is the eurrent state, x^^^ is the next state, and is the ehange in states from 
one iteration to the next. The ehange in states is described as 

^k=A^,-A=KPk (A.2) 

where the vector is the search direction and the scalar \ is the step size. Equation 
(A.2) ean be rearranged to for iteration to give the next state as 

A^i=^k+\Pk (A.3) 

In order for seareh algorithm to eonverge, the next value of the performanee function 
must be less than the eurrent for eaeh iteration, sueh as 

^ )<yiA) for \ (A.4) 

with the seareh direction vector, , being less than zero and the step size being small, 
but greater than zero [85]. 

There are various ways for seleeting both the seareh direetion and the step size. 
The Steepest Deseent method searehes in the direction of the largest negative rate of 
change for the performance funetion. The step size for each iteration ean be either fixed 
or variable. The performanee funetion is assumed to be an analytie function, such that its 
derivatives exist. This is logieal sinee our performanee functions are based on position 
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states, so the first derivative is veloeity, and the seeond derivative is aeceleration. Any 
performance function could be represented by a Taylor Series expansion, such as 


V{x) = V{x*) + —V{x) 
dx 


(x-x*) + 






dx‘ 


-V{x) 


{x-x*f +. 


(A.5) 


where the state x is a scalar and the function is expanded about some nominal point, x *. 
For the neighborhood of the nominal point, the higher order terms can be dropped for a 
quadratic estimate of the function about that point. Since the states are most likely a 
vector, such as a position vector, the Taylor series can be expanded into a vector notation 
with y(T) =y(Xj,X 2 ,...x„), as 


y(T)=y(T*)+—y(T) 

ox. 


d 


(x,-x,*)+—y(T) 


dx^ 


(^2 -^2 *) + ••• 


+—V(T) 

ox„ 


1 


d^ 


1 /-\ 



U; 



(y-^i*)" + 


(A.6) 



d^ 



\2) 

0^2 


x=x^ 

f 1 A 


(X 2 -X 2 *)"+...+ 






-y(T) 




Since this format is cumbersome, it is more convenient to write in matrix form 




y(T) = y(T )+vy(T) (t-t )+ — (t-t ) vy (T) 

\x=x* 

where Vy is the gradient of the potential function, and is defined as 

vy(T) = 


(A.7) 


-^V{x),^V{x),...,^V{x) 


dx. 


dx. 


dx„ 


(A.8) 


and Vy^ is the Hessian, and is defined as 
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VV\x) = 


5" 


y(^) 


dx^dx2 


V(x). 


dx^dx^ 


V(x) 


5^2 5^1 


V(x) 


dXr, 


-V{x)... 


dx^dx^ 


-V{x) 


v{x) -—y(^)... ^v(^) 


dx^dx^ 


Sx„dx2 


dx„ 


(A.9) 


The Hessian is the Jaeobian of the gradient, where the Jaeobian of a function, f(x), is 
defined as a matrix of its first partial derivatives 


Jix) = 


dxj 

A 

A 


fl(x) 

A 

A 

flix) 

A 

A 


ox^ 



(A.10) 


The gradient and the Hessian can give valuable information about the performance 
surface, which may also be referred to as the potential field. The gradient is the direction 
of maximum slope from the current position. On the performance surface, any direction 
which is orthogonal to the gradient will have zero slope [85]. 

The baseline quadratic performance function is of the form 

V{x) = ^{xY Ax + d^x + c (A.ll) 


where x are the states of the function, A is a scaling matrix which relates to the shape of 
the contour ellipses, J is a scaling factor which determines where equilibrium point is 
located, and c is a scaling constant which changes the value of the performance function 
at the equilibrium point. The gradient and Hessian of a generic quadratic function are 

VV{x) = Ax + d (A. 12) 

VV\x) = A (A.13) 
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The equilibrium point is generally seleeted to be at the origin (d = 0) with a minimum 
value equal to zero (c = 0), so the quadratie funetion simplifies as 

Vix) = ^ixyAx (A.14) 

In this case, the gradient and Hessian of the simple quadratic become 

VV{x) = Ax (A.15) 

VV\x) = A (A. 16) 

Notice that both the gradient and the Hessian depend on the performance surface contour 
shaping matrix, A. 

(A.17) 

Using this information about the gradient and substituting into the general performance 
search algorithm in Equation (A.3), yields the Steepest Descent algorithm. 

(A. 18) 

where the search vector, , is in the direction of the negative gradient, -g ^,, which 
relates as 

A=-vy(x)|^_, (A.19) 

IX —Xk 

The Steepest Descent algorithm for the simple quadratic becomes 

(A.20) 

The Steepest Descent algorithm applied to a simple quadratic gives insight into 
the optimization schemes characteristics. For a simple quadratic performance surface, 
with A being an identity matrix and the state values being the relative distance from the 
equilibrium point, the performance surface is shown in Figure A.l. In this simple case 
the minimum is located at the origin and the contour lines of equal potential are perfect 
circles. An iterative search along the steepest gradient can be shown to converge, but the 
convergence path will vary dependent on the step size parameter. 
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S-axis R-axis 

Figure A.l Sample Performanee Surfaee. 


Determination of the step size parameter, \ ,ean be eonstant or variable. 

Constant step size parameters are general based on a desired convergenee response. 
Evaluating the largest eigenvalues of the Hessian ean ensure eonvergenee and even 
provide an underdamped response. For example, consider a system described by the 
quadratic function with 


2 0 
A = 

0 50 


(A.21) 


and its minimum at the origin. This function is steeper along the y-axis then the x-axis 
and has eigenvalues of (2, 50). The maximum eigenvalue, of the quadratic 

function can be used to determine the bounds on the step size parameter. For 
convergence, the step size must meet the following condition. 


0 < 2 -< 



(A.22) 
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For an overdamped response, the step size must meet the following eondition, 


0 < 1 < 


r 1 ^ 
V^^'^max J 


(A.23) 


The eonvergenee eharaeteristies of four different step sizes, using the funetion 
deseribed in Equation (A.21), are shown with the 3D performanee surfaee in Figure A.2 
and as eontours in Figure A.3. For the step size meeting the eonditions of Equation 
(A.23), /I = 0.01, the Steepest Descent method converges with an overdamped behavior. 
For the step size meeting the condition of Equation (A.23), /I = 0.03, the Steepest 
Descent method converges with some overshoot behavior. For the case where the step 
size is equal to the upper range condition of Equation (A.23), X = 0.04, the method does 
not converge and continues to oscillate. If the step size is increase beyond this condition, 
X = 0.041, the method diverges. Notice that the search directions are always opposite the 
gradient and thus perpendicular to the contour lines. 
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Figure A.2 Steepest Descent Performance for Fixed Step Sizes in 3D. 
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Figure A.3 Steepest Descent Contour for Fixed Step Sizes. 


Variable step sizes may be based on minimization-on-a-line or momentum 
learning methods [85]. Minimization-on-a-line allows variation of the step size at each 
iteration, k . The analytical minimization for a quadratic function is satisfied by setting 
its derivative equal to zero and solving for the step size [85]. 


\=- 


f - r - A 
SkPk 

Pk A p, 


(A.24) 


The Hessian for simple quadratic functions is symmetric, so the transpose may be 
dropped in these cases. The convergence characteristics of the variable step size, using 
the function described in Equation (A.24), are shown with the 3D performance surface in 
Figure A.4 and as contours in Figure A.5. The convergence generally takes less 
iterations, but dramatically oscillates due to the changing of the search direction at each 
iteration. Each step is orthogonal to the previous step and parallel to every other step. 
This pattern is due to the primary axis of the contour shape and the tangential intercept of 
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the contour lines. This direction change may not be reasonable with systems which have 
large momentum and can not change direction rapidly. 



Figure A.4 Steepest Descent Performance for Variable Step Sizes in 3D. 
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Figure A.5 Steepest Deseent Contour for Variable Step Sizes. 


Momentum learning smoothes out the iterations and may speed eonvergenee for 
eonsistent trajeetories [85]. This is done by eonsidering both the eurrent and previous 
iteration information and weighting them to reduee oseillations. The momentum learning 
eoeffieient, /ul , aeting as a low-pass filter on the seareh direetion oseillation, weighs the 
new seareh direetion with respeet to the previous seareh.. The Steepest Deseent with 
fixed step size and momentum learning yields the following algorithm 

[(1 - aO gk - (aO gk-i] (A.25) 

For jul = 0, the funetion has no memory of previous searehes and is the original Steepest 
Deseent algorithm. For jul = 1, the seareh direetion is not updated from iteration to 
iteration. Momentum learning is not generally useful for variable step size 
implementation. The step size is determined from seareh direetions whieh are damped by 
the momentum and as a result overshoot usually oeeurs. This may inerease the number 
of iterations neeessary for eonvergenee. For those interested, a diseussion of both 
Newton’s Method and Conjugate Gradient seareh algorithms follows. 
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B. NEWTON’S METHOD 

The Newton’s Method is based on loeal seeond order approximation, as an 
extension of the Steepest Deseent first order approximation. The Newton’s Method 
algorithm is 

-Sk (A.26) 

For simple quadratic functions, Newton’s method converges in one step. This is 
dependent on the Hessian being positive definite. The Newton’s method usually 
converges faster than the Steepest Descent method. However, the second order 
approximation is computationally expensive, due to the need to calculate a Hessian and 
its inverse. 

A first order variation, referred to as the Gauss-Newton Method, only requires the 
calculation of the Jacobian matrix, J , and its inverse. The Gauss-Newton algorithm is 

^^+1 ) (A.27) 

where v(T^) comes from rewriting Equation (A.l 1) as V(x) = v{xf *v{x). This method 
neglects second order terms and can result in the Hessian not being invertible. 

Another variation, called the Levenberg-Marquardt Scheme, converges faster at 
the cost of additional memory. This is allows for a single parameter adjustment, juk , to 
vary the search performance from a Gauss-Newton scheme to the Steepest Descent with a 
small step size. The Levenberg-Marquardt algorithm is 

+ (A.28) 

where juk is determined so that the inverse matrix exists. As juk ^0 the algorithm 
approaches the Gauss-Newton algorithm and as juk the algorithm approaches the 
Steepest Descent with a small learning rate. This method converges fast, but has a 
drawback of requiring more memory. 

C. CONJUGATE GRADIENT 

Since using a second order approximation may be too expensive, a first order 
approach called conjugate gradient may be effective. This approach steps along the 
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conjugate directions of the function, with the initial step being in the direction of negative 
gradient. The step size can be adjusted to minimize the function along the selected search 
direction. For larger problems the conjugate gradient method is more efficient, since it 
does not require computation of the eigenvalues of the Hessian or its inverse. It usually 
converges in a number of interactions which is equal to the number of states. For 
instance a simple three state, [x, y, z], quadratic function would converge in three steps. 
However, as the function varies from quadratic then the iteration scheme will require 
more steps to converge. A primary conjugate gradient algorithm of interest is called the 
Fletcher Reeves iteration scheme. This iteration scheme iterates as follows 


^k^i=^k+\-A 

(A.29) 


(A.30) 

Pk^i=-Sk^i+A^i-Pk 

(A.31) 

f ^ ^ ^ ^ 


P Sk+i Sk+i 

^k+i ~ -r- 

(A.32) 

1 j 



The step size, \, varies as minimization along a line, refer to Equation (A.24). 
The gradient, g ^., is calculated as previously discussed. However, now the search 
direction, , varies from the Steepest Descent after the first iteration. The scalar 
variation of search direction, 5^, can take may forms, but all are based on Gram-Schmidt 
orthogonalization concepts [85]. Notice that the initial search direction is along the 
Steepest Descent, = -gg, such that 5^ = 0 . 

The convergence characteristics of the conjugate gradient Fletcher Reeves 
iteration scheme, using the function described in Equations (A.29) - (A.32), are shown 
with the 3D performance surface in Eigure A.6 and as contours in Eigure A.7. The search 
direction of the two iteration steps is slightly greater than 90 degrees. This illustrates the 
variation of the search direction which allows the conjugate gradient approach to 
converge quicker then the Steepest Descent method, even though they are using the same 
variable step size. 
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Figure A.6 Conjugate Gradient with Fleteher Reeves Seheme. 
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Figure A.7 Conjugate Gradient Contour with Fleteher Reeves Seheme. 
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D. DISCUSSION 

The Steepest Deseent seareh method is the simplest and most straight forward of 
the diseussed methods. The Newton’s Method and Conjugate Gradient approaehes 
require additional computation. In the dynamics environment the additional computation 
may not yield any improvement in the control algorithm performance. However, if the 
control algorithm cycle is at least twice as fast as the measurement cycle then 
implementation of the Conjugate Gradient search algorithm should be considered. 
Practical limitations in actuator performance may show that the desired changes in search 
directions are impractical. For instance, drastic directional changes are impossible for 
physical system with momentum and limited actuation, such as spacecraft. Although, the 
same might be observed in a Steepest Descent algorithm which does not decrease in step 
size as it approaches equilibrium point. 

Modification may also be made in the state weighting to influence the initial and 
final search directions. This could be useful for constrained docking approaches where 
convergence along a relative axis is necessary. The weighting of the states should have 
direct correlation with the cost and performance functions discussed previously. By 
modifying the contour shape of the quadratic surface the search algorithm can be 
prejudiced to approach from the desired direction. These types of modification usually 
negatively influence efficiency due to the forcing of required, but artificial, maneuver. 
This has a direct comparison with the collision avoidance capability. 

It is not envisioned that modifications to the search algorithm will dramatically 
improve the APF convergence performance. Although, some improvement may be 
possible in control efficiency or maneuver duration. It was initially thought that 
replacing the state scaling function, in Equation (A.30), with the Clohessy-Wiltshire’s A 
matrix from Equation (4.16) may allow for the linearized spacecraft dynamics to be 
incorporated into the APE. This relationship appears to have merit when used for the 
spatial states, but requires further analysis. In the APE development velocity 
considerations were the primary driving states for both successful convergence and 
collision avoidance. Determination of desirable velocity parameters may become much 
less intuitive for various close proximity maneuvers. So, simple reliance on the 

geometrically based APE may be limited regardless of the search method applied. 
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APPENDIX B: LQR WITH COLLISION AVOIDANCE GAINS 


During the course of developing the LQR/APF with collision avoidance control 
algorithm, consideration was given to the direct incorporation of collision avoidance 
parameters into the LQR gain matrixes. It was believed to be a logic step to consider the 
inclusion of collision avoidance weighting into the gain matrixes of the LQR’s cost 
function; refer to Equation (7.1). Modification to the gain matrixes for the states, Q ,and 
control effort, R , continued to be of primary interest. There was some initial success in 
selecting additive functions that resulted in successful collisions avoidance, however 
rigorous evaluation is still necessary. For instance, it is unclear what effect additive gain 
terms modification may have on stability, convergence, and efficiency. 

The approach was to add positive semi-definite parameters into the state gain 
matrix without increase the number of states. Although adding relative position and 
velocity states for each obstacle could conceptually work, including these full or partial 
states for each potential obstacle causes the computational requirements to increase 
dramatically. To avoid this consequence, the inclusion of obstacle information might be 
accomplished by increasing the cost of core spacecraft’s states in the established relative 
structure. This should influence the control effort to avoid obstacle regions. 

Numerous attempts to include additive terms in the gain matrixes were made. 
Additive positive, semi-definite terms were only considered, since the gains matrixes 
must maintain the proper form. It was determined that relative position information 
could not be incorporated with only additive terms, due to the complications of weighting 
position terms of obstacles in the presence of the goal. However, the inclusion of 
additive velocity terms, , showed potential. The additive velocity term is based on the 

Gaussian shaping function, from Equation (7.19), and the magnitude of the Chase 

spacecraft’s velocity component in the direction of the obstacle, |v^^| based on Equation 

(7.17). If the Chase spacecraft velocity is positive in the relative direction of the obstacle 
then the additive velocity gain terms was determined to be 
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where the small positive constant, , is included to prevent numerical difficulty due to 
division by zero. The additive velocity gain magnitude shaping function, was 

selected to increase influence in response to velocity toward obstacles while decreasing in 
influence in the vicinity of the goal. This decrease in influence is a scaled version of the 
safety function, k^, determined in Equation (7.21). The resulting magnitude shaping 
function is 
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The additive terms for the velocity gain components are included along the diagonal of 
the LQR state weighting matrix and modify Equation (7.7) as follows 


Q = 


— 00 
r 

eg 

0 — 0 


0 

0 

0 


r 

eg 


0 0 — 

r 

eg 


0 

0 

0 


0 


0 

0 

0 


0 

i^ + CJk^ 
0 
0 


0 


0 

(1 + C'v2)^e 


0 


0 


0 

0 

(1 + C.3)^e 


(B.3) 


where k„ was determined from Equations (7.10) and (7.12) to be 
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No modification was made to the eontrol effort weighting matrix, R , sinee eontrol effort 
should not be limited in the presence of obstaeles. The ineorporation of terms that would 
decrease weighting was not considered due to the challenges of maintaining the proper 
eost relationship between various states. 

The LQR additive eollision avoidanee gains were generally sueeessful in avoiding 
obstaeles and eonverging toward the goal position. Preliminary results show that the 
LQR with eollision avoidanee gains eontrol algorithm could avoid obstacles and other 
spaeecraft while eonverging to within 4.0 meters of the Target spaceeraft’s outer 
boundary. The results for the six close proximity maneuvers with collision avoidance are 
listed in Table B.l. For these eollision avoidanee maneuvers, the goal position is the 
eenter of the Target spaceeraft. This requires that the Target spaeeeraft’s repulsion to 
allow the Chase spaeeeraft to eonverge while avoiding impaet. Stationary obstaeles are 
plaeed at positions along the unobstrueted path of the Chase spaeeeraft. These obstaeles 
have an actual diameter of 2.0 meters. Generally, avoiding larger obstacles requires more 
eontrol effort and time. The performanee aceuraey shows that general convergence and 
rendezvous ean be aehieved while enabling collision avoidance. The relative position, 
veloeity and eontrol effort plots of the first Chase spaeecraft are shown in Figure B.l 
through Figure B.3. Although densely paeked obstaele regions tend to have an additive 
repulsion foree, whieh keep later eonverging spaeeeraft at further distanees away from 
the mutual goal loeation. For instanee, the three Chase spaeeeraft in the near rendezvous 
maneuvered to within 1.0, 2.0, and 4.0 meters of the Target spaeeeraft, based on their 
order of arrival. 


Rendezvous Maneuver 

LQR with CA Gains 

Near with Obstaele 

Av = 0.6760 m/s 

Initial RSW [0, 70, 0] m 

td= 1210 s 

Near with Obstacle 

Av = 0.5196 m/s 

Initial RSW [50,-100, -50] m 

td= 1073 s 

Near with Obstaele 

Av = 0.6845 m/s 

Initial RSW [100, 100, 100] m 

td= 1069 s 


Table B.l Rendezvous Maneuver with Collision Avoidance Results. 
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Velocity Difference in RSW Coordinates, m/s 



Figure B. 1 Chase Spaeeeraft Relative Position. 



Figure B.2 Chase Spaeeeraft Relative Velocity. 
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There were several diffieulties and defieiencies this initial approaeh. First, Chase 
spaeecraft approaching along the relative coordinate frame’s W-axis did not properly 
avoid obstacles directly along their path. In these cases the Chase spacecraft slowed in 
the region of the obstacle, but continued to push toward converges. If the obstacle gain 
was increased. Chase would avoid collision by stopping and no longer converging. This 
problem may preventable by perturbing the Chase off of this perfectly balanced 
approach. This could be done by adding a shifting factor if an obstacle is detected 
directly along this W-axis path. Also, adjusting the weighting the W-axis motion could 
influence the Chase spacecraft to move along the orbital plane. 

Next, the focus on velocity terms may lead to some chattering when a Chase 
spacecraft skirts along an obstacle’s region of influence. In this situation the relative 
velocity toward the obstacle fluctuates on the exterior of the avoidance arch. These 
fluctuations in close proximity would result in large cost variation. This causes a 
chattering phenomenon, as shown in Figure B.3. This issue may require additional gain 
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terms or control logic to be included. The velocity gains as initially set to increase near 
obstacles, but may need more complex functional relationships as they evade in the 
region of influence. 

Finally, docking in the presence of obstacle was not successful. Problems 
developed when state cost decrease during goal convergence could not be properly tuned 
with collision avoidance in this terminal region. If too large, the collision avoidance 
terms prevented convergence; if too small, the obstacles were not properly avoided. The 
superposition of numerous obstacles in a region also, resulted in challenges for this 
conceptual approach. These problems were previously overcome by applying decaying 
obstacle influence along the geometric region of approach. The application the APF 
based techniques did not work for this approach. The cost variations would change the 
feedback response in a fashion that would restrict the usefulness of the logic. For 
instance, if a the obstacle’s gain parameters increase the cost function to much, such that 
the Chase spacecraft is not longer approaching, then damping of that function only slows 
divergence. This leads to local minimums around some obstacles. The obstacle 
avoidance parameters may be shaped to avoid these cost circumstances away from the 
goal position, but there is less possible cost variation in the region of the Target 
spacecraft. 

The high degree accuracy needed for more precision rendezvous and docking 
maneuvers was not achieved, so direct comparison with Chapter VI results could not be 
made. Attempts to balance the gains of position and velocity in the vicinity of the Target 
spacecraft tended to result in chattering. Additional research may still prove successful in 
the incorporation of additive collision avoidance terms in the LQR’s state gain matrix. 
Any success will need to take complex weighing and state relationships into 
consideration. These relationships may prove more complicated then practically 
realizable. The benefits of such an algorithm are not yet apparent. 
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